)OT/FAA/SE-91/4 

)OT-VNTSC-FAA-91-13 

IAS  System 
ingineering  Service 
i/ashington,  DC  20591 


AD-A242  323 


Multisensor  Signal 
Processing  Techniques 

(Hybrid  GPS/LORAN-C  with  RAiM) 


Dr.  Frank  van  Graas 
Mark  R.  Kuhl 

Ohio  University 
Athens,  OH  45701 


Final  Report 
September  1991 


This  document  is  available  to  the  public 
through  the  National  Technical  Information  Service 
Springfield,  Virginia  22161 


91-14316 


U.S.  Department  of  Transportation 

Federal  Aviation  Administration 

9110  2« 


NOTICE 


This  document  is  disseminated  under  the  sponsorship  of  the 
Department  of  Transportation  in  the  interest  of  information 
exchange.  The  United  States  Government  assumes  no 
liability  for  its  contents  or  use  thereof. 


NOTICE 

The  United  States  Government  does  not  endorse 
products  or  manufacturers.  Trade  or  manufacturers' 
names  appear  herein  solely  because  they  are  considered 
essential  to  the  object  of  this  report 


1.  Raport  No.^  \J 

DOT/FAA/SE-91/4 

2.  Gevarnment  Accession  No. 

4.  Titl#  and  Subtitl# 

Multisensor  Signal  Processing  Techniques 
(Hybrid  GPS/LORAN-C  with  RAIM) 

Df.^l^anlc  van  Grass,  Mark  R. 

Kuhl* 

9.  Parlorfntng  Organisotion  Nomo  ond  Address 

Avionics  Engineering  Center 
Ohio  University 

Athens,  OH  45701 

Tvchnieol  Report  Documontotion  Page 


3.  R«cipi«nf’s  Cotoiog  No. 


5-  Rupert  Oav* 

September  1991 


Organisation  Cod* 


6.  Porlerming  Organisotion  Roport  No. 

DOT-VNTSC-FAA-91-13 


10.  Work  Unit  No.  (TRAIS) 

FA160/A1Q74 _ 

It.  Contpoct  OP  Gpont  No. 

DTRS-57-87-C-000006 


13.  Typo  of  Ropopt  ond  Popiod  Covopod 

Final  Report 

August  1989  -  October  1990 


14.  SponioPing  Agoney  Codo 

ASE-300 


12.  Spentoping  Agoney  Nomo  ond  Addpoti 

U.S.  Department  of  Transportation 
Federal  Aviation  Administration 
NAS  System  Engineering  Service 
Washington,  DC  20591 


15.  Suppiomontopy  Notoi  U  .  ii .  Department  oT^ransportation 

*Under  Contract  to;  Research  and  Special  Programs  Administration 

John  A.  Volpe  National  Transportation  Systems  Center 
Cambridge,  MA  02142 


16.  Abttpoet 

One  of  the  major  elements  in  alleviating  existing  problems  in  en  route  airspace  is 
to  allow  more  aircraft  to  traverse  a  given  volume  of  airspace.  Recent  developments 
in  navigation  systems  will  support  this  effort  by  enabling  user  preferred  routes  and 
by  offering  more  precise  and  reliable  navigational  capabilities. 

Multisenor  navigation  system  architectures  are  presented  with  the  emphasis  on  signal 
processing  techniques  and  receiver  autonomous  integrity  monitoring  (RAIM) .  This  is 
followed  by  a  detailed  study  of  a  hybrid  navigation  receiver  based  on  the  NAVSTAR 
Global  Positioning  System  (GPS)  and  the  Long  Range  Navigation  System  LORAN-C.  The 
design,  implementation  and  flight  testing  of  a  prototype  hybrid  GPS/LORAN  receiver  are 
presented,  including  the  realtime  simulation  and  detection  of  navigation  signal 
malfunctions . 


17.  K*y  Wspdf 

Radionavigation  System,  Radio 
Determination  System,  Position 
Location,  Global  Positioning 
System,  Long  Range  Navigation, 
Fault  Detection 


19.  Saewhty  (o1  this  roport) 

UNCLASSIFIED 


18.  Oiltributien  Statsmant 


DOCUMENT  IS  AVAILABLE  TO  THE  PUBLIC 
THROUGH  THE  NATIONAL  TECHNICAL 
INFORMATION  SERVICE.  SFHINGFIELO. 
VIRGINIA  33161 


20.  Saeurity  Claatif.  (of  ihis  pago) 

UNCLASSIFIED 


21.  No.  of  ragoa  22.  Fnca 

96 


Form  DOT  F  1700.7  (8-72) 


Ropreduetion  of  complatad  poga  outliarisad 


PREFACE 


The  work  desscribed  in  this  final  report  was  performed  in  support  of  the 
Federal  Aviation  Administration.  Jerry  W.  Bradley  (ASE-300)  is  the  sponsor's 
program  manager. 

This  report  is  a  deliverable  under  Task  2  of  PPA  No.  FA- 160.  The  work 
was  performed  by  Assistant  Professor  Frank  van  Graas  of  Ohio  University, 
Athens,  Ohio.  The  PPA  manager  is  M.  J.  Moroney  (RSPA/VNTSC/DTS-52) . 

One  of  the  major  elements  in  alleviating  existing  problems  in  en  route 
airspace  is  to  allow  more  aircraft  to  traverse  a  given  volume  of  airspace. 
Recent  developments  in  navigation  systems  will  support  this  effort  by  enabling 
user  preferred  routes  and  by  offering  more  precise  and  reliable  navigational 
capabilities.  One  develpment  of  specific  interest  is  the  integration  of  GPS 
and  LORAN-C  to  acheive  a  sole  means  of  navigation.  This  report  addresses  the 
capabilities  offered  by  multisensor  navigation  systems  in  general  with 
emphasis  on  signal  processing  techniques  and  receiver  autonomous  integrity 
monitoring  (RAIM) .  In  addition  this  report  details  the  design, 
implementation,  and  flight  testing  of  a  protype  hybrid  GPS/LORAN  navigation 
system. 
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One  of  the  major  elements  in  alleviating  existing  problems  in  en  route 
airspace  is  to  allow  more  aircraft  to  traverse  a  given  volume  of  airspace. 
Recent  developments  in  navigation  systems  will  support  this  effort  by  enabling 
user  preferred  routes  and  by  offering  more  precise  and  reliable  navigational 
capabilities.  One  development  of  specific  interest  is  the  integration  of  GPS 
and  LORAN-C  to  achieve  a  sole  means  of  navigation.  This  report  addresses  the 
capabilities  offered  by  multisensor  navigation  systems  in  general  with 
emphasis  on  signal  processing  techniques  and  receiver  autonomous  integrity 
monitoring  (RAIM) .  In  addicion  this  report  details  the  design,  implementation, 
and  flight  testing  of  a  protype  hybrid  GPS/LORAN  navigation  system. 

The  first  part  of  this  report  shows  the  derivation  of  the  ordinary  j east 
squares  (OLS)  estimator  from  the  linear  model  using  the  Projection  Theorem. 

The  performance  of  the  OLS  is  subjected  to  a  detailed  error  analysis  which 
includes  the  effects  of  geometry  and  measurement  errors,  both  bias  and  noise. 
Using  similar  procedures,  the  extended  Ridge  and  Kalman  filters  are  derived 
which  are  representative  for  the  classes  of  biased  and  unbiased  estimators, 
respectively.  It  is  shown  that  in  the  presence  of  measurement  bias  errors, 
both  estimators  converge  to  the  OLS  inflation  of  the  position  error.  It  is 
therefore  concluded  that  in  the  case  of  dominant  measurement  bias  errors  (e.g. 
GPS  or  LORAN) ,  the  integrity  monitoring  performance  offered  by  the  OLS 
estimator  cannot  be  improved  upon.  It  was  also  found  that  the  Ridge  estimator 
can  be  used  to  explain/optimize  the  performance  of  a  mismatched  Kalman  filter. 
This  provides  a  very  helpful  insight  into  the  performance  of  the  Kalman  filter 
in  the  presence  of  unmodeled  dynamics. 

Based  on  the  conclusion  that  the  OLS  estimator  is  sufficient  to  perform 
integrity  monitoring,  a  general  solution  methodology  is  presented  for  a 
multisensor  navigation  system.  The  range  residual  in  egrity  parameter  was 
derived,  again  using  the  Projection  Theorem. 

The  above  results  have  been  applied  to  a  realtime  prototype  hybrid 
GPS/LORAN  receiver.  Two  flight  tests  revealed  that  the  hybrid  GPS/LORAN 
receiver  performs  in  accordance  with  its  design.  The  course  deviation 
indicator  is  responsive  and  the  indicated  course  compares  favorably  with  other 
area  navigation  equipment.  Out  of  a  total  of  fifteen  simulated  signal 
malfunctions,  twelve  malfunctions  which  would  have  caused  unacceptable  course 
deviations  were  detected  by  the  range  residual  integrity  algorithm. 

Equal  weighting  of  GPS  and  LORAN  measurements  is  used  for  the  flight 
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tests.  Therefore,  the  accuracy  of  the  hybrid  system  will  be  mostly  determined 
by  the  LORAN  measu^'ements .  Standard  LORAN  propagation  models  are  used  such 
that  the  achieved  accuracies  are  represer .ative  for  current  LORAN  receivers. 
Because  of  this,  the  accuracy  of  the  hybrid  system  will  not  be  as  good  as  that 
provided  by  GPS;  however,  the  availability  and  integrity  of  the  hybrid  system 
exceeds  that  of  GPS  by  several  orders  of  magnitude  [4].  At  the  same  time,  the 
hybrid  system  accuracies  were  still  found  to  b?  well  within  all  current 
requirements  as  listed  in  Chapter  3.  It  is  also  concluded  that  the  accuracy 
of  hybrid  GPS/LORAN  can  be  improved  upon  significantly  through  LORAN 
calibration  or  by  using  a  we.ighting  matrix  in  the  hybrid  solution,  see  Seccion 
6.1. 


In  summary,  hybrid  GPS/LORAN  is  a  successful  case  study  of  fully  hybrid, 
multisensor  navigation.  Similar  performance  characteristics  may  be 
anticipated  for  other  multisensor  systems  based  on  sensors  such  as  GPS.  LORAN, 
GLONASS,  Omega,  and  baro-altimeter . 
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INTRODUCTION 


1 . 


One  of  the  major  elements  in  alleviating  existing  problems  in  en  route 
airspace  is  to  allow  more  aircraft  to  traverse  a  given  vol\ime  of  airspace  [1]. 
Recent  developments  in  navigation  systems  will  support  this  effort  by  enabling 
user  preferred  routes  and  by  offering  more  precise  and  reliable  navigational 
capabilities.  Included  in  these  developments  are  new  satellite  technologies 
such  as  the  NAVSTAR  Global  Positioning  System  (GPS)  and  improvements  of 
existing  systems  such  as  LORAN-C.  One  development  of  specific  interest  is  the 
integration  of  GPS  and  LORAN-C  to  achieve  a  sole  means  of  navigation.  Hybrid 
GPS/LORAN  is  setting  the  stage  for  a  next-generation  of  multisensor  navigation 
systems  which  is  anticipated  to  have  a  major  impact  on  the  safety, 
reliability,  and  efficiency  of  national  and  invernational  transportation.  A 
highly  reliable  multisensor  navigation  system  will  increase  low-altitude, 
remote  and  offshore  operations.  Additionally,  a  multisensor  navigation  system 
could  serve  as  the  position  sensor  for  automatic  dependent  surveillance  (ADS) 
and  Communication,  Navigation,  and  Surveillance  (CNS)  services,  thereby 
increasing  the  safety  and  operational  efficiency  of  the  Air  Traffic  Control 
(ATC)  System. 

The  first  part  of  this  report,  Chapters  2  through  6,  addresses 
multisensor  navigation  systems  in  general,  with  emphasis  on  signal  processing 
techniques  and  receiver  autonomous  integrity  monitoring  (RAIM) .  RAIM  is 
recommended  to  be  used  as  the  integrity  technique  for  a  multisensor  navigation 
system.  The  second  part  of  this  report.  Chapters  7  through  10,  details  the 
design,  implementation,  and  flight  testing  of  a  prototype  hybrid  GPS/LORAN 
navigation  system.  Chapters  11  and  12  contain  conclusions  and  recommendations 
for  future  efforts,  respectively. 
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2. 


PRINCIPLES  OF  POSITION  FIXING 


2.1  Path  Delay  Systems 

Two  basic  methods  exist  for  a  vehicle  to  determine  its  position:  dead 
reckoning  and  position  fixing.  Position  determination  using  a  dead  reckoning 
method  involves  projecting  a  known  position  forward  to  some  future  time  by 
keeping  continuous  track  of  velocity,  accelerations,  or  both.  Examples  are 
Doppler  and  Inertial  Reference  Systems.  A  position  fix  method  involves  a 
vehicle  sensing  landmarks  and  determining  an  instantaneous  position  without 
reference  to  any  former  position.  Landmarks  are  at  known  locations  on  or 
around  the  earth.  Examples  of  landmarks  are  the  Very  High  Frequency 
Omnidirectional  Range  (VOR)  stations,  the  Global  Positioning  System  (GPS) 
satellites,  and  stars  (celestial  fixes). 

A  vehicle  can  obtain  range  (p)  or  bearing  (6)  information  from  the 
landmarks  by  measuring  parameters  such  as  phase  delay  and  pulse  time  of 
arrival  (TOA)  delay.  These  make  up  what  are  called  path  delay  systems  (see 
Kelly  and  Cusick  [2]).  Extending  the  classifications  presented  by  Kelly  and 
Cusick,  path  delay  systems  may  be  of  four  different  types:  single  path, 
multiple  paths  in  three  dimensions,  multiple  paths  in  four  dimensions,  and 
combinations.  The  four  different  types  of  path  delay  systems  are  shown  in 
Figure  2.1. 

Single  path  systems  may  be  classified  as  either  active  or  passive.  An 
active  single  path  system  requires  only  one  landmark  to  determine  a  range. 

The  user  in  an  active  single  path  system  transmits  a  signal  to  the  landmark 
and  listens  for  a  reply.  The  range  is  then  determined  by  multiplying  half  the 
round  trip  delay  time  by  the  speed  of  transmission  of  the  signal.  To  measure 
the  round  trip  delay  time,  the  user  clock  needs  only  to  be  stable  during  the 
measurement  time.  A  crystal  oscillator  is  a  sufficient  frequency  reference 
source  for  an  active  single  path  measurement.  Examples  of  active  single  path 
systems  include  the  Distance  Measuring  Equipment  (DME)  system  and  radio 
detecting  and  ranging  (RADAR)  systems.  A  receiver  in  a  passive  single  path 
system  measures  path  delay  times  (or  ranges)  with  respect  to  known 
transmissions  from  landmarks.  This  method  requires  that  the  user  clock  be 
synchronized  with  the  transmitter.  Furthermore,  both  the  user  clock  and  the 
transmitter  clock  must  be  stable  throughout  the  mission  time.  As  an  example, 
if  one  of  the  frequency  reference  sources  drifts  one  part  in  10^^  (Rubidium 
standard),  then  after  1000  seconds  the  clock  would  have  drifted  10 
nanoseconds,  which  corresponds  to  a  drift  in  the  range  measurement  of 
approximately  10  feet.  Passive  single  path  systems  are  currently  not  used  for 
vehicle  navigation. 

Three-dimensional  multiple  path  systems  measure  differences  in  reception 
time  (or  compare  the  path  delays)  between  at  least  two  landmarks  and  the 
receiver,  or  use  multiple  range  measurements.  These  systems  may  be  classified 
as  true  hyperbolic,  radial  or  spherical.  True  hyperbolic  systems  determine 
when  time  or  frequency  differences  in  the  path  delays  are  constant.  This 
yields  a  hyperbolic  line  of  position  (LOP).  When  the  time  or  frequency 
difference  in  the  path  delays  are  zero,  a  straight  LOP  will  result.  This  is 
called  a  radial  system.  Examples  of  hyperbolic  systems  are  LORAN-C  and  Omega, 
while  an  example  of  a  radial  system  is  VOR.  An  example  of  an  active  spherical 
system  is  multiple  DME.  Spherical  systems  reduce  to  a  circular  system  if  the 
altitude  of  the  vehicle  is  known.  An  example  of  a  passive  two-path  system  is 
the  two-station  French  LORAN-C  chain  (3]. 
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Four-dimensional  multiple  path  systems  are  similar  to  spherical  systems 
without  synchronized  clocks.  The  passive  measurements  to  the  landmarks  are 
called  "pseudoranges"  (or  "pseudofrequencies")  because  the  measurements 
contain  an  unknown  clock  phase  (or  frequency)  offset.  Therefore,  the  user 
must  solve  for  its  clock  phase  or  frequency  offset  with  respect  to  the 
landmarks.  The  unknown  receiver  clock  phase  or  frequency  adds  another 
dimension  to  be  solved  for  in  the  position  solution.  Examples  of  four¬ 
dimensional  multiple  path  systems  are  the  Global  Positioning  System  (GPS)  and 
the  Global  Navigation  Satellite  System  (GLONASS).  These  two  worldwide 
satellite  systems  are  sponsored  by  the  United  States  and  the  Soviet  Union, 
respectively. 

In  order  to  improve  accuracy,  availability,  reliability  and  integrity, 
the  above  path  delay  systems  may  be  combined  [4] .  With  the  focus  on  aircraft 
positioning,  combinations  of  the  various  path  delay  systems  can  be  implemented 
in  two  ways ; 

(1)  Multi-Sensor  System 

(2)  Hybrid  System 

A  multi-sensor  system  incorporates  navigation  data  from  landmarks  which 
are  selected  by  the  system  operator  or  by  some  type  of  automatic  selection 
scheme.  As  specified  in  Federal  Aviation  Administration  (FAA)  Advisory 
Circular  20-130  [5],  a  multi-sensor  system  may  be  approved  for  instrument 
flight  rules  (IFR)  enroute  and  terminal  navigation  within  the  conterminous 
United  States  (CONUS).  But  for  instrument  approaches,  only  navigation  data 
from  landmarks  which  are  already  individually  approved  for  charted  instrument 
approaches  may  be  selected.  Moreover,  a  sole  means  of  navigation  must  be 
installed  and  operable  at  the  same  time  [6], 

A  hybrid  system  uses  navigation  data  from  all  landmarks  that  are 
available  for  all  phases  of  flight  (enroute,  terminal  and  approach).  For  a 
hybrid  system,  no  selection  process  is  exercised.  An  example  of  a  hybrid 
system  may  be  GPS/LORAN-C.  Alone,  GPS  and  LORAN-C  will  not  meet  sole  means 
requirements  for  navigation  [7].  But  by  combining  these  two  systems  into  a 
hybrid  system,  sole  means  requirements  could  be  satisfied  [8] .  A 
"multisensor"  hybrid  system  -  or  simply  called  a  multisensor  system  -  arises 
from  the  fact  that  many  systems  may  be  combined  for  navigation,  i.e.  a  hybrid 
system  consisting  of  GPS,  GLONASS,  DME,  VOR,  Omega  and  LORAN-C.  Currently, 
navigation  requirements  have  not  been  defined  for  hybrid  systems,  see  Chapter 
3. 


2.2  Classification  of  Position  Fix  Methods 

The  way  a  vehicle  determines  a  position  fix  from  the  information 
provided  by  the  landmarks  may  be  classified  into  several  different  methods: 

(1)  Triangulation  Method 

(2)  Trilateration  Method 

(3)  Polar  Method 

(4)  Hyperbolic  Method 

(5)  Combinations  of  Methods 

These  methods  are  discussed  below  for  a  two-dimensional  (2-D)  position  fix. 
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(1)  Triangulation  Method  (6-6) 


Bearing  information  from  two  different  landmarks  may  be  used  to  obtain  a 
position  fix.  This  method  is  called  triangulation  and  is  depicted  in  Figure 
2.2.  The  method  proceeds  as  follows;  Each  of  the  two  stations  provides  the 
vehicle  with  a  certain  bearing  angle  from  true  or  magnetic  North  with  respect 
to  the  station.  These  bearing  angles  will  cross  at  only  one  location,  giving 
a  unique  position  fix.  VOR  stations  are  an  example  of  a  navigation  system 
that  transmits  bearing  information  to  the  vehicle. 

(2)  Trilateration  Method  (p-p) 

Range  information  from  two  different  landmarks  may  be  used  to  obtain  a 
position  fix.  This  method  is  called  trilateration  and  is  portrayed  in  Figure 
2.2.  The  method  proceeds  as  follows;  A  particular  range  is  provided  by  each 
of  the  two  stations.  Each  range  defines  a  circle  of  positions  on  which  a 
vehicle  may  be  located.  The  two  circles  will  cross  at  two  locations  -  one 
being  the  correct  vehicle  position  and  the  other  being  an  ambiguous  solution. 
The  ambiguous  solution  may  be  eliminated  by  using  prior  information, 
additional  information  such  as  speed  or  heading,  or  a  third  station.  DME 
stations  are  an  example  of  a  navigation  system  that  provides  range  information 
to  the  vehicle.  In  three  dimensions,  the  range  from  a  DME  station  is  called  a 
slant  range  to  take  into  account  the  altitude  of  the  vehicle.  Also,  GPS  and 
GLONASS  are  navigation  systems  from  which  a  vehicle  can  determine  a  three- 
dimensional  position  fix  based  on  ranges  between  the  vehicle  and  four 
satellites  in  the  system.  (The  fourth  satellite  is  needed  to  resolve  the 
receiver  clock  offset  from  the  system  clock.) 

(3)  Polar  Method  (p-e) 

Range  and  bearing  information  from  a  single  landmark  may  also  be  used  to 
acquire  a  position  fix.  This  method  is  called  the  polar  method  and  is  shown 
in  Figure  2.3.  The  method  proceeds  as  follows;  One  station  provides  the 
vehicle  with  both  a  bearing  angle  and  a  range.  As  stated  before,  the  range 
gives  a  circle  of  possible  positions;  and  the  bearing  angle  indicates  where  on 
that  circle  the  vehicle  is  located.  The  VOR/DME  stations  and  Tactical  Air 
Navigation  (TACAN)  stations  are  examples  of  navigation  systems  that  transmit 
both  range  and  bearing  information  to  the  vehicle. 

(4)  Hyperbolic  Method 

Another  method,  called  the  hyperbolic  method,  uses  lines  of  constant 
time  difference  (TD)  between  two  pairs  of  landmarks  to  determine  a  position 
fix.  The  hyperbolic  method  is  also  depicted  in  Figure  2.3;  and  it  requires 
three  different  landmarks  (stations  1  and  2  make  up  one  pair  while  stations  2 
and  3  constitute  the  second  pair).  The  method  proceeds  as  follows;  Station  1 
transmits  a  signal  which  is  received  by  the  user  at  some  later  time.  The  same 
transmission  is  also  received  by  station  2.  Upon  reception  of  the  signal  from 
station  1,  station  2  transmits  a  signal.  A  constant  difference  between  the 
time  of  arrivals  of  the  signals  from  station  2  and  station  1  define  a 
hyperbolic  LOP,  on  which  the  user  is  located.  In  the  same  way,  a  second  pair 
of  stations  is  used  to  determine  another  hyperbolic  LOP.  Where  these  two 
hyperbolic  LOPs  cross  gives  a  position  fix  for  the  vehicle.  LORAN-C  and  Omega 
are  examples  of  navigation  systems  that  employ  time  differencing  for  position 
fixing. 
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STATION  1 


STATION  2 


,ure  2.2  Triangulation  and  Trilateration  position  fix  methods 


As  a  special  note,  LORAN-C  may  be  used  in  the  passive  ranging  mode  as 
well  as  in  the  hyperbolic  mode  of  operation. 


(5)  Combinations  of  Methods 

The  above  methods  may  be  combined  in  various  ways  to  obtain  a  position 
fix.  This  gives  rise  to  a  multisensor  navigation  system:  a  vehicle  may 
integrate  information  from  all  the  available  navigation  systems  to  obtain  an 
optimum  position  fix. 


2,3  Position  Fixing  Accuracies  of  Current  Path  Delay  Systems 

Range  and  bearing  measurements  are  subject  to  a  variety  of  error 
sources,  including  measurement  errors,  noise  sources  (atmospheric  and  man¬ 
made),  multipath,  unmodeled  ionospheric  and  tropospheric  propagation  delays, 
clock  errors  and  landmark  location  uncertainty.  (References  [9]  and  [10] 
provide  a  detailed  description  of  the  error  sources  affecting  the  GPS  in 
particular. ) 

Measurement  errors  propagate  into  position  errors  as  a  function  of 
geometry  and  solution  methods.  Table  2.1  summarizes  the  accuracies  of  current 
U.S.  radionavigation  systems  (compiled  from  the  1990  edition  of  the  Federal 
Radionavigation  Plan  (FRP)  [11]).  These  accuracies  are  95X  confidence  level 
numbers,  meaning  that  statistically  at  least  95X  of  the  position  estimates 
will  be  within  the  value  listed.  An  exception  to  this  accuracy  level 
specification  are  the  2-D  accuracies.  The  2  drms  (distance  root  mean  square) 
numbers  are  used  where  2-D  accuracies  are  specified.  As  shown  in  reference 
[12],  2  drms  numbers  may  have  anywhere  from  a  952  to  a  982  confidence  level. 
(Two  different  accuracies  are  shown  for  the  GPS  in  Table  2.1  because  the  GPS 
offers  two  navigation  services:  Precise  Positioning  Service  (PPS)  for 
military  users  and  Standard  Positioning  Service  (SPS)  which  anybody  can  use.) 
At  this  time,  the  accuracies  for  GLONASS  are  not  clearly  defined,  although  the 
system  should  have  similar  performance  characteristics  as  the  GPS. 
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Navigation  System 


Accuracy  (95Z  confidence  level) 


VOR 

±  1.4° 

DME 

185  meters 

TACAN 

±  1°,  185  meters 

Radiobeacons 

±3-10° 

Omega 

3.7  -  7.4  kilometers 

LORAN-C 

460  meters 

Transit 

single  frequency 

500  meters 

dual  frequency 

25  meters 

GPS 

SPS  (civil) 

100  meters  (horizontal) 
156  meters  (vertical) 

PPS  (military) 

17.8  meters  (horizontal) 
27.7  meters  (vertical) 

Table  2.1  Accuracies  of  current  navigation  systems  as 
specified  by  the  1990  FRP  [H]- 
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3. 


REQUIREMENTS  FOR  AIRBORNE  POSITION  FIXING  SYSTEMS 


Two  types  of  navigation  systems  must  be  considered  when  addressing  the 
navigation  requirements:  1.  sole  means  of  navigation  systems,  and  2. 
supplemental  type  navigation  systems.  The  main  difference  between  these  two 
types  of  systems  is  that  the  sole  means  of  navigation  system  must  have  some 
level  of  redundancy  to  allow  for  continued  navigation  in  the  presence  of 
signal  failures.  Supplemental  systems  only  have  to  announce  the  failure,  upon 
which  the  pilot  switches  to  a  sole  means  of  navigation.  Therefore, 
differences  in  system  requirements  are  mainly  in  the  areas  of  integrity  and 
availability. 

Five  major  performance  characteristics  that  must  be  satisfied  are 
accuracy,  availability,  reliability,  coverage  and  integrity  (see  reference 
[11,  13,  lA]  for  additional  requirements).  Furthermore,  the  system  must  be 
compatible  with  the  air  traffic  control  system  and  the  existing  airspace 
structure,  which  is  largely  based  on  VOR/DME.  Most  of  the  requirements  for 
sole  means  navigation  are  well  understood  and  dociunented,  except  for 
availability  and  integrity  requirements.  The  FAA  is  currently  developing  a 
regulation  that  will  establish  minimum  standards  under  which  a  radionavigation 
system  may  be  certified  as  the  sole  radionavigation  system.  The  Advance 
Notice  of  Proposed  rulemaking  has  been  distributed  [13]. 

The  focus  of  this  report  is  on  the  feasibility  of  receiver  autonomous 
integrity  monitoring  schemes  to  supply  the  integrity  information.  In  the  case 
of  a  multisensor  navigation  system,  many  redundant  measurements  are  available 
to  determine  the  validity  of  the  position  solution. 

For  supplemental  navigation,  the  performance  requirements  in  the  areas 
of  accuracy  and  integrity  are  taken  from  reference  [15],  and  are  summarized  in 
Table  3.1.  The  following  notes  are  provided  to  explain  some  of  the  quantities 
in  Table  3.1: 

1.  Detection  probability  is  defined  as  a  conditional  probability;  it  is 
the  probability  of  not  detecting  an  out-of-tolerance  condition  given 
that  an  out-of-tolerance  condition  is  present.  The  detection 
probability  must  be  guaranteed  at  each  space/time  point.  It  is  assumed 
that  errors  that  are  hard  to  detect  are  slowly  growing  bias  errors  (of 
order  2  m/s)  and  that  these  errors  occur  approximately  once  every  two 
years  for  each  satellite,  which  corresponds  to  a  probability  of 
occurrence  of  10*®.  This  then  yields  an  overall  detection  probability 
of  1  -  10*®.  Consequently,  the  alarm  threshold  would  be  exceeded  with  a 
probability  of  10*’.  (The  detection  probability  is  on  a  per  sample 
basis ) . 

2.  All  accuracy  specifications  are  on  a  95Z  probability  basis. 
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Phase  of  Flight 


Enroute 


Terminal 


Nonprecision 

Approach 


Performance  |  111 

- 1 - 1 - 1 - 1 

- T 

Position  Fixing  1 

Radial  Accuracy  | 

1 

0.12A  nmi  1 

1 

- T 

0.124  nmi  | 

1 

1 

0.124  nmi 

- 1 

1 

1 

1 

Flight  Technical  | 

Error  (crosstrack)  | 

1 

1.00  nmi  1 

(2.00  oceanic)  j 

1 

1 

1.00  nmi  j 

1 

1 

0.50  nmi 

1 

1 

1 

Total  Radial  | 

Accuracy  j 

1 

1 

1.008  nmi  | 

(2.004  oceanic)! 

1 

1.008  nmi  | 

1 

.  J. 

0.515  nmi 

1 

1 

-  I  ■ 

Radial  Alarm  Threshold] 

I 

 -  j-. 

2e00  nmi  | 

1 

- h 

1 .00  nmi  j 

1 

0.30  nmi 

=1 

1 

_ 1 

I 

Maximum  Allowable  ( 

Alarm  Rate  | 

...  1 

0.002 /Hr  i 

1 

1 

0.002 /Hr  i 

1 

1 

0.002 /Hr 

1 

1 

1 

Time  to  Alarm  | 

1 

1 

30  sec*  1 

1 

10  sec.  1 

! 

10  sec. 

1 

1 

t 

Minimum  Detection  | 

Probability  j 

- L 

.999  i 

1 

_ L 

r 

.999  1 

1 

1 

.999 

1 

1 

_j 

Table  3.1  Tentative  accuracy  and  integrity  requirements 
for  supplemental  GPS  (from  ref.  [15]). 
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4. 


THEORY  OF  LINEAR  REGRESSION 


This  section  covers  the  general  topic  of  estimating  a  vehicle  position 
fix  given  range  measurements  with  respect  to  landmarks.  The  ordinary  least 
squares  (OLS)  estimator  is  presented;  and  it  is  shown  that  the  OLS  estimator 
results  in  an  unbiased  position  solution.  An  error  criterion,  called  the  mean 
square  error  (MSE),  is  defined  in  order  to  analyze  the  performance  of  position 
estimators.  Next,  the  effects  of  geometry  on  the  position  fix  are  addressed. 
Finally,  a  DME  simulation  is  presented  which  provides  insight  into  the 
performance  of  the  OLS  estimator  under  poor  geometry  conditions. 


4.1  The  General  Linear  Model 

As  stated  in  Section  2.1,  all  path  delay  systems  provide  either  a  range 
or  bearing  from  which  a  vehicle  can  obtain  a  position  fix.  The  position  fix 
is  an  estimate  because  of  errors  in  the  measurement  data.  Since  the  data 
exhibit  a  significant  degree  of  error  or  "noise",  the  strategy  is  to  find  a 
position  fix  that  best  represents  the  general  trend  of  the  data,  which  is 
accomplished  through  the  use  of  regression  theory  [16]. 

For  this  section,  assume  that  the  measurement  data  consist  of  range 
measurements,  which  are  related  to  a  position  fix  by  a  nonlinear  model.  Since 
the  relation  is  nonlinear,  one  must  first  linearize  the  relation  by  taking  a 
Taylor  series  expansion  around  some  given  reference  point  (usually  the  initial 
estimate)  [17],  see  Section  4.4.1.  The  general  linear  model  is  given  by 

Y=Hfi+AB+e  (4.1) 

where  Y  is  the  n  x  1  data  vector  which  consists  of  n  range  measurements,  H  is 
the  known  n  x  p  design  matrix  which  may  include  linearized  terms  which 
represent  a  nonlinear  system  and  is  the  p  x  1  regressor  vector  to  be 
estimated.  For  navigation,  the  elements  of  g  may  be  position,  velocity, 
acceleration  and  clock  offset,  (g  is  also  known  as  the  system  state  vector.) 
The  e  vector  represents  zero  mean,  uncorrelated,  normally  distributed  noise  on 
the  range  measurements,  i.e.  e  -  N(0,R)  where  R  is  the  known  covariance  matrix 
of  e  and  N  signifies  the  normal  distribution.  The  covariance  matrix  is 
determined  as  follows 

R  -  E[ee'^]  -  0^1  (4.2) 

where  E  signifies  the  expected  value,  a  is  the  standard  deviation  of  the 
measurement  noise  and  I  is  the  identity  matrix  (the  off-diagonal  elements  of  R 
are  zero  because  the  noise  is  uncorrelated).  Bias  errors,  represented  by  the 
AB  vector,  also  exist  on  the  range  measurements.  Bias  errors  are  seen  as 
constant  and  unknown  to  the  estimator.  In  this  section,  the  analysis  of 
noise-only  on  the  measurement  data  will  be  performed.  Section  4.5  covers  the 
case  of  both  noise  and  bias  errors.  Without  bias  errors,  the  linear  model  is 
given  by 


Y  -  Hg  +  e  (4.3) 

Many  techniques  have  been  developed  to  solve  the  above  equation;  the  most 
straightforward  being  the  ordinary  least  squares  (OLS)  estimator. 
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U .2  Ordinary  Least  Squares  (OLS)  Estimator 

The  OLS  estimator  is  developed  as  follows.  The  residuals  after 
estimating  in  the  linear  model  are  given  by 

Sr  =  Y  -  Hi  =  Y  -  Y  (4.4) 

where  e^  is  the  residual  vector,  i  is  an  estimate  of  g  and  Y  is  the  estimated 
measurement  vector.  Using  the  Projection  Theorem  [18,19],  project  the 
measurements,  Y,  from  the  measurement  space  onto  the  estimation  space  such 
that  the  distance  between  true  Y  and  estimated  Y  is  a  minimum.  This  is  shown 
in  Figure  4.1.  Therefore,  Y  is  orthogonal  to  ejj,  or 

-  0 

(Hi)T(Y  -  Hi)  -  0 
i'^H^CY  -  Hi)  -  0 

i'^(H'^Y  -  H'^Hi)  «  0  (4.5) 

The  problem  is  to  solve  this  equation  for  i.  Since  i  is  usually  not  zero,  the 
terms  inside  the  parenthesis  must  be  set  to  zero.  This  yields  the  OLS 
estimator 


H^Y  -  H^Hi  -  0 

H^Y  =  H'^tii  -*  called  the  "normal  equation" 


i  =  (hTh)-ihTy  =  ioLs  (^-6) 

Substituting  equation  4.3  into  equation  4.6  gives 
ioLs  “  (HTH)-lHTHfi  +  (H^H)-lH^e 

lots  “  ^  (^-7) 

From  equation  4,7  the  expected  value,  or  mean,  of  the  OLS  estimator  is 

E[ioLsJ  =  E[fi  +  (H^H)*lH’'e] 

^  (^*8) 

since  the  expected  value  of  e  is  zero.  Therefore,  the  OLS  estimator  is  an 
unbiased  estimator.  The  bias  is  zero,  denoted  by 

BIAS[fioLsl  -  0  (^-9) 

The  variance  of  the  OLS  estimator  is 

VAR[1olsJ  =  E[goLS^]  -  [E(goLs)l^ 

VAR[ioLsJ  “  o2tRACE[(H^H)-^)  (4.10) 


where  a  is  the  standard  deviation  of  the  measurement  noise  and  the  TRACE 
operator  is  the  sum  of  the  diagonal  elements  of  a  matrix. 


13 


H,  H, 


^'lODEL  ESTIMATE 


TRUE  MODEL 


I  i 

T  T 

\  h  -1  h  ■; 

I  ll  *5]  ll 


MEASUREMENT  VECTOR 


TRUE  MODEL 
MEASUREMENT  NOISE 


OLS  RESIDUALS 


ESTIMATE  ERROR 


C2X  G29  /  005  E58 


Figure  4.1  Measurement  space  and  estimation  space  with  Projection 
Theorem  geometry  (reference  [19]). 


4.3  Mean  Square  Error  (MSE)  Error  Criterion 

The  OLS  estimator  minimizes  the  sum  of  t..  .  squares  of  the  residuals, 
SS(ep),  which  is  given  by 

SS{eR)  =  (Y  -  Hl)’’(Y  -  Hi)  (4.11) 

where  i  =  ioLS  minimum  SS(ej^).  However,  for  navigation  purposes  the 
appropriate  error  criterion  is  the  mean  square  error  (MSE).  The  MSE  indicates 
how  far  off  the  estimated  position  is  from  the  true  position.  It  includes 
both  bias  errors  in  the  solution  as  well  as  variance.  For  the  OLS  estimator, 
the  MSE  is  given  by 

MSE[ioLsl  =  EtdoLS  -  -  fi)] 

MSE[ioLs]  “  +  (HTH)-iHTe  -  ^)T(g  +  (HTH)-lHTe  -  £)] 

MSE[ioLs]  “  OUTRAGE [(H'^H)-M  (4.12) 

which  is  the  same  as  the  variance  (given  by  equation  4.10).  This  should  be 
expected  since  the  OLS  estimator  is  an  unbiased  estimator,  no  bias  exists  in 
the  solution.  Another  way  to  determine  the  MSE  is  as  follows 

MSE[^QLg]  =  BIAS[^ql5]^  +  VAR[^QLg] 

MSE[goLs]  -  0  +  VAR[ioLs] 

MSE[^ols1  “  OUTRAGE [(H'^H)-M  (4.13) 

In  the  case  of  vehicle  positioning,  the  OLS  estimator  minimizes  the 
SS(eD),  not  the  MSE.  Since  the  OLS  estimator  does  not  allow  a  bias  to  exist 
in  the  solution,  it  restricts  the  MSE  to  include  only  the  variance  term.  This 
variance  inflates  significantly  when  a  "poor  geometry"  condition  arises  among 
the  landmarks  and  the  vehicle.  Geometry  effects  on  position  fixing  will  be 
discussed  in  the  next  section. 


4.3.1  Geometry  Effects  on  Position  Fixing 

For  position  fixing,  the  geometry  of  the  landmarks  with  respect  to  the 
vehicle  is  a  major  factor  in  how  accurate  the  position  estimate  will  be  in 
terms  of  the  MSE.  In  order  to  illustrate  what  is  meant  by  the  term  geometry, 
equation  4.12  is  converted  into  its  canonical  coordinates  using  a  singular 
value  decomposition  (SVD)  of  H 

gnxp  .  ynxn  j.nxp  (4.14) 

where  U  and  are  defined  as  orthogonal  transformation  matrices  that 
diagonalize  H,  and  2  is  a  diagonal  matrix  whose  elements  are  the  singular 
values  of  H.  The  singular  values  of  H  are  equal  to  the  square  root  of  the 

eigenvalues  of  H  (vX|).  For  a  more  in  depth  discussion  of  the  SVD  and  how  it 
relates  to  solving  OLS  problems,  see  reference  (20J .  Therefore,  equation  4.12 
transforms  into 

MSEf^oLS^  "  O^TRAGEIS'^Z] 
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(A. 15) 


^^1  J 

where  are  the  eigenvalues  of  h’^H.  Note  that  as  the  eigenvalues  become 
smaller,  the  MSE  grows.  Small  eigenvalues  represent  an  ill-conditioned  (or 
nearly  singular)  H  matrix.  This  gives  rise  to  what  is  considered  to  be  a  poor 
geometry  condition. 

Figure  4.2  shows  a  physical  interpretation  of  the  relation  between  the 
eigenvalues  and  the  position  fix  in  two  dimensions.  When  estimating  a  two- 
dimensional  position  (2-D),  a  statistical  error  ellipse  is  defined  [21].  (the 
error  ellipse  will  contain  95Z  to  98Z  of  the  position  estimates,  see  reference 
[12].)  This  error  ellipse  is  a  function  of  the  range  measurement  noise  and 
the  eigenvalues.  Assuming  that  the  measurement  noise  is  fairly  constant,  the 
error  ellipse  will  inflate  when  the  eigenvalues  become  small.  For  the  2-D 
case,  one  of  the  two  eigenvalues  will  obtain  a  small  value  as  the  crossing 
angle,  7,  becomes  small. 

For  example,  y  “  90 '  is  considered  to  be  very  good  geometry  while  y  = 

1°  is  considered  to  be  very  poor  geometry.  The  crossing  angle  is  also 
depicted  in  Figure  4.2.  If  the  range  measurements  from  the  two  DME  stations 
in  Figure  4.2  have  a  small  crossing  angle,  X2  will  be  small;  and  this  will 
inflate  the  error  ellipse  in  the  direction  corresponding  to  X2. 

"Near  collinearity"  is  another  term  that  is  used  for  the  case  when  a 
small  crossing  angle  exists  between  the  range  measurements  [22]  .  Nearly 
collinear  range  measurements  inflate  the  range  errors  into  a  poor  position 
estimate.  This  is  what  is  meant  by  a  poor  geometry  condition. 

4.3.2  GDOP  Factor  Description 

For  navigation,  the  type  of  geometry  condition  that  exists  is  portrayed 
by  a  statistical  factor  called  the  Geometric  Dilution  of  Precision  (GDOP). 

The  GDOP  is  defined  as  follows 


GDOP  -  VTRACE(H'^H)*^  (4.16) 

Five  different  types  of  DOPs  are  defined  below 

GDOP  -  3-D  position  and  time  dilution  of  precision 
PDOP  “  3-D  position  dilution  of  precision 

HDOP  ■  position  dilution  of  precision  in  the  horizontal  plane 
VDOP  -  position  dilution  of  precision  in  the  vertical  direction 
TDOP  =  time  dilution  of  precision 

For  satellite  navigation  (GPS,  GLONASS),  GDOP  is  used  to  represent  a  4- 
D  geometry  condition:  3-D  position  and  time  (X,Y,Z,T).  In  the  context  of 
navigation  using  the  GPS,  Jorgenson  [23]  states  the  following  about  GDOP 

1)  GDOP  is,  in  effect,  the  amplification  factor  of  pseudorange 
measurement  errors  into  user  errors  due  to  the  effect  of 
satellite  geometry. 

2)  GDOP  is  independent  of  the  coordinate  system  employed. 
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3)  GDOP  is  a  criterion  for  designing  satellite  constellations. 

4)  GDOP  is  a  means  for  user  selection  of  the  four  best  satellites 
from  those  that  are  visible. 

PDOP  is  more  useful  in  navigation  than  GDOP  because  it  is  the  critical 
information  needed  to  calculate  position.  Also,  one  may  wish  to  calculate 
PDOP  by  finding  the  volume  of  a  tetrahedron  rather  than  implementing  equation 
4.16  [24].  As  seen  in  reference  [24],  this  saves  processing  time  over  the 
matrix  inversion  method  of  equation  4.16  when  many  possible  combinations  of 
satellites  exist. 

Furthermore,  for  aircraft,  the  vertical  position  information  is  obtained 
from  the  baro-altimeter .  The  good,  relative  accuracy  of  the  baro-altimeter  is 
the  basis  for  the  vertical  separation  of  aircraft.  Therefore,  aircraft 
positioning  using  landmarks  only  requires  2-D  information,  and  subsequently, 
HDOP  is  the  parameter  of  interest.  Note  that  HDOP  is  dependent  on  the 
coordinate  system  used  as  it  involves  the  propagation  of  measurement  errors 
into  a  user-defined  horizontal  plane. 


4.4  DME  Simulation  Example 

A  DME  simulation  is  presented  to  illustrate  the  points  made  above  about 
geometry  conditions  and  how  they  affect  the  OLS  estimator.  The  simulation 
involves  the  estimation  of  an  aircraft’s  position  as  it  travels  along  a 
constant  velocity  flight  path.  Figure  4.3  shows  the  simulation  set-up. 


4.4.1  DME  Simulation  Set-Up 

Range  measurements  from  two  different  DME  stations  are  needed  to  solve 
for  a  two-dimensional  position  (X,Y)  as  depicted  in  Figure  4.3.  The  equations 
are  shown  below 

[(X-Xi)2  +  (Y-Yi)2]l/2  -  (4.17) 

where  Rj  is  the  range  measurement  from  the  aircraft  to  DME  station  i,  (X,Y)  is 
the  estimated  user  position  and  (Xj^,Yj^)  are  the  coordinates  for  DMEj^.  Using 
these  two  equations,  the  aircraft  estimates  its  position  and  velocity.  One 
way  to  do  this  is  to  linearize  the  above  equation  by  constructing  a  Taylor 
series  expansion  around  an  initial  estimate  and  retaining  only  the  first  order 
terms . 


aRj 

+  - 

6X  + 

— 

6Y 

(4.18) 

ax 

(Xq.Yo) 

aY 

(Xo.Yq) 

where  6X  and  6Y  are  the  corrections  to  the  initial  aircraft  state  estimate 
(Xq,Yq).  R^q  is  the  initial  estimated  range  to  DME  station  1  using  the 
initial  aircraft  state  estimate.  This  is  shown  below 


[(Xo-Xi)2  +  (Yo.Yi)2]1/2  -  R^o 


(4.19) 
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The  direction  cosines •  ^iji  cosines  of  the  angles  between  the 

line-of-sight  (LOS)  vector  from  the  aircraft  to  DME  station  i  and  its 
projection  along  the  coordinate  axes 


X  -  X. 


Y  -  Y. 


(A. 20) 


The  measurement  model  for  a  constant  velocity  flight  path  is  given  below 
for  a  batch  of  j  measurement  sets;  j  »  1,  N.  N  is  the  number  of  range 
measurements  in  a  batch  interval  (AT^).  See  also  Figure  4.3. 
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(4.21) 


where  i  is  the  DME  station  number,  j  is  the  measurement  set,  and  Vy  is  the 
aircraft  velocity,  and  6Rij  ■  R^j  -  Rio*  is  the  intersample  time  of  the 

range  measurements  in  a  batch. 

6X  and  6Y  do  not  represent  the  difference  in  position  from  one 
observation  to  the  next.  They  represent  the  initial  aircraft  position  offset 
from  the  start  of  the  constant  velocity  flight  path;  or  in  other  words,  the 
error  in  the  a  priori  estimate  of  the  aircraft  position  (Xq,Yq).  The 
direction  cosines  are  calculated  every  batch-interval  using  new  estimates  of  X 
and  Y  (X,Y). 

Based  on  one  batch  of  measurements,  the  aircraft  estimated  state  is 
updated  at  time  tj^  as  follows 


+  (N-l)AT, 

6Y  V, 


(4.22) 


where  Vy  and  Vy  are  the  estimates  of  velocity. 


Results  of  the  OLS  Estimator  Using  the  DME  Simulation 

In  order  to  evaluate  the  performance  of  the  batch  OLS  estimator  for 
navigation,  two  cases  are  anal7zed: 

(1)  A  good  geometry  situation:  7  «  90  (HDOP  “  1.4). 

DME  station  1  location,  [Xj  Yj)  “  [-l.lOSeOS  0]  feet 

DME  station  2  location,  [Xj  Y2]  “  (3.608e05  0]  feet 

(2)  A  poor  geometry  situation;  7-1°  (HDOP  -  80). 

DME  station  1  location,  [Xj  Yj]  “  [3.464e05  0]  feet 

DME  station  2  location,  [Xj  Y2j  “  [3.608e05  0)  feet 

Both  cases  have  the  following  initial  conditions; 

true  aircraft  location,  [X^j^  Y^j^]  «  [0  2.000e05)  feet 

measurement  noise,  o  ■  60  feet 

measurement  bias,  ^  -  [00]  feet 

number  of  range  measurements  in  a  batch,  N  -  16 

aircraft  initial  position  offset,  [6X  6Y]  •  (270  -370]  feet 

For  batch  estimation  with  measurement  noise  only,  note  that  the  position  error 
inflation  is  effectively  reduced  by  one  over  the  square  root  of  the  number  of 
measurements  in  a  batch.  Since  a  batch  of  16  measurements  is  used  in  the 
simulation,  the  position  error  inflation  above  is  reduced  by  a  factor  of  4. 

For  these  two  cases,  a  1000  iteration  Monte  Carlo  analysis  is  performed. 
For  this  simulation,  the  first  iteration  of  the  simulation  process  is  repeated 
M  times  (M  *  1000),  and  each  time  different  measurement  noise  is  used.  The 
results  for  cases  (1)  and  (2)  are  shown  in  Figures  4.4  and  4.5,  respectively. 
Notice  that  the  errors  in  Figure  4.4  are  in  a  near  circular  formation.  This 
is  due  to  the  fact  that  good  geometry  exists  (7  ■  90°).  Figure  4.5  depicts 
the  position  results  of  the  OLS  estimator  when  a  poor  geometry  condition 
exists  (7  “  1°);  the  solution  is  still  unbiased,  but  the  errors  form  a  rather 
large  ellipse  around  the  true  solution.  The  MSE,  which  is  the  appropriate 
error  criterion  to  try  to  minimize  for  navigation  purposes,  is  given  below  in 
terms  of  the  square  root  of  the  MSE  (i.e.  SQRT[MSE])  for  each  case: 

(1)  SQRT(MSE]  -  20  feet  (2)  SQRT[MSE]  -  1177  feet 

From  these  results  it  can  be  seen  that  the  OLS  estimator’s  MSE  and 
associated  error  ellipse  is  inflated  significantly  when  a  poor  geometry 
situation  exists.  To  see  just  when  the  MSE  starts  to  inflate  as  the  geometry 
becomes  worse,  the  results  from  "small"  Monte  Carlo  analyses  for  various 
crossing  angles  between  the  measurements  are  compiled.  In  order  to  obtain 
different  crossing  angles,  the  location  of  DME  station  1  is  moved  in  distinct 
steps  from  the  good  geometry  location  ([-1.108e05  0]  feet)  toward  the  DME 

station  2  location  ([3.608e05  0]  feet).  Each  step  is  4716  feet  in  the  X 

direction.  At  each  step  a  50  iteration  Monte  Carlo  analysis  is  implemented 
and  the  SQRT[MSE]  is  recorded.  The  crossing  angle  versus  the  SQRT[MSE]  is 
shown  in  Figure  4.6.  From  this  figure  it  can  seen  that  the  SQRT[MSE]  starts 
to  increase  rapidly  when  the  crossing  angle  becomes  less  than  6  . 
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Fig’jre  4,4  Monte  Carlo  analysis  of  the  OLS  estimator  in  a  good 
geometry  condition  (7  =  90°). 


22 


^  -lOOOh 

■o 


-2000 


-3000 L 


-4000  ' - ' - - - - - ■ — 

-3000  -2000  -1000  0  1000 

5X  Error  (feet) 


2000 


Figure  4.5  Monte  Carlo  analysis  of  the  OLS  estimator  in  a  poor 

geometry  condition  (-y  =  1°). 


4.5  OLS  Estimator  with  Range  Measurement  Bias  Errors 

To  complete  the  discussion  on  the  OLS  estimator,  the  MSE  equation  must 
be  developed  for  the  case  of  bias  errors  as  well  as  noise  existing  on  the 
range  measurements.  Start  with  the  linear  model  given  by  equation  4.1 

Y-H^+AB+e  (4.23) 

Recall  that  the  OLS  estimator  is  given  by 

Sols  “  (H^H)-^H^Y  (4.24) 

Substituting  the  linear  model  into  the  OLS  estimator  yields 

Sols  *  (H’^H)  *  ^H'^HS  +  (H'^H)  '  ^H'^AB  +  (H'^H)-^H'^e 

loLS  '  ^  +  (HTH)'^H^AB  +  (H’‘H)-^H^e  (4.25) 

From  equation  4.23  the  expected  value,  or  mean,  of  the  OLS  estimator  is 

E[fioLs3  “  E[S  +  (H^H)'*h'^A6  +  (H’‘H)-^H^e] 

E[fioLsl  “  ^  +  (H^H)-^h’’M  (4.26) 

Therefore,  the  OLS  estimator’s  solution  bias  is 

BIAS(Sols1  “  (H''‘H)'^H'^AB  (4.27) 

This  shows  that  the  bias  of  the  OLS  solution  is  only  a  function  of  the 
measurrjent  bias  errors.  If  no  measurement  bias  errors  exist,  then  the  OLS 
solution  would  be  unbiased,  as  was  shown  before  in  Section  4.2.  The  variance 
of  the  OLS  estimator  is  the  same  as  before 

VAR[£ols]  -  -  CE(fioLs)3^ 

VAR[ioLsJ  =  a2TRACE[(H%)-l]  (4.28) 

The  MSE  equation  is  then  given  by 

MSE[fioLs]  -  EK^ols  -  -  S)] 

MSECfioLsJ  "■  E[(fi  +  (H^H)-*H’^AB  +  (H’'H)'^H^e  -  £)'^ 

(g  +  (H’''H)*^h’'AB  +  (H^H)-^H^e  -  fi)  J 

MSE[ioLsJ  ”  Ab'''H(H'^H)-^'‘'(H''’H)-^h'^AB  +  OUTRAGE  [  (H'^'H)  '  ^ 

Since  (H^H)'^  is  symmetric,  (H^H)'^^  *  (H^H)'^  and  the  MSE  becomes 

MSEEioLsl  “  AB^H(H‘^H)*2h''‘AB  +  OUTRAGE [  (H^'H)  '  h  (4.29) 

This  completes  the  discussion  on  the  OLS  estimator,  which  will  be  used  as  the 
baseline  estimation  scheme  to  evaluate  the  performance  of  the  Kalman  filter 
and  the  Ridge  regression  signal  processor  presented  in  the  next  Section. 
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5. 


THE  EXTENDED  KALMAN  AND  RIDGE /KALMAN  FILTERS 


5.1  Theory  of  Kalman  Filtering 

This  section  moves  from  the  previous  discussions  on  batch  processing  to 
a  recursive  processing  technique.  The  Kalman  filter  is  presented,  which  can 
be  shown  to  be  just  a  recursive  form  of  the  ordinary  least  squares  (OLS) 
estimator.  For  two  key  references,  [25]  portrays  the  inherent  equivalence  of 
OLS  estimation  and  Kalman  filter  theory;  and  [26]  shows  that  Kalman  filter 
theory  is  essentially  the  same  as  random  parameter  linear  regression  theory. 


5.1.1  Recursive  versus  Batch  Processing 

The  performance  of  the  batch  and  the  recursive  processor  is  equivalent. 
Batch  processing  gives  a  position  estimate  after  averaging  a  group  of 
measurements  over  a  specified  time  interval.  Usually,  it  is  considered 
memory less,  since  it  only  analyzes  the  measurement  group  over  the  given  time 
interval;  no  weight  is  given  to  the  measurements  that  were  received  before  the 
time  interval  began.  A  batch  processor  may  also  incorporate  some  type  of 
previous  information  (called  a  priori  information),  weighted  appropriately, 
into  the  linear  model  [27]  .  On  the  other  hand,  a  recursive  processor  employs 
a  time-history  of  measurements  received  as  the  a  priori  information,  and 
continuously  updates  the  current  position  estimate  (or  state)  with  a 
combination  of  the  newly  received  measurement  data  and  the  model  estimate 
obtained  from  the  previous  measurement  information. 

Below  it  will  be  shown  how  the  Kalman  filter  is  developed  from  the  model 
with  a  priori  information.  Next,  issues  are  discussed  involving  the  selection 
of  the  error  covariance  matrix  (Pq)  and  system  covariance  matrix  (Q) . 


5.1.2  Kalman  Filter  Development 

In  1960,  R.  E.  Kalman  provided  an  alternative  way  of  formulating  the  OLS 
estimator  using  state-space  methods  [28].  The  result  is  what  is  known  today 
as  the  Kalman  filter.  There  are  many  ways  to  derive  the  Kalman  filter 
including;  the  Bayesian  approach,  using  the  Projection  Theorem,  and  from 
properties  of  Gaussian  conditional  probability  densities.  Below,  the  Bayesian 
approach  is  implemented  to  develop  the  Kalman  filter  from  a  model  that 
includes  a  priori  information.  This  will  be  shown  to  form  the  Kalman  filter’s 
updating  step.  Since  the  Kalman  filter  is  a  recursive  estimator,  its  updating 
step  will  allow  new  information  to  be  combined  in  some  way  with  previously 
given  measurement  information  at  each  time  interval. 

Assume  that  the  linear  model  given  by  equation  4.1  in  Section  4.1  has 
experimental  data  concerning  such  that  the  average  of  all  previous 
measurements  equals  jSq  with  some  uncertainty,  v.  Using  and  v,  two 
equations  now  define  the  model;  the  initial  conditions  equation  and  the 
linear  model  equation 

S.  -  e.Q  +  V  (5.1) 

X  -  +  AB  +  e  (5.2) 

Recall  that  ^  represents  range  measurement  bias  errors  and  e  -  N(0,R)  defines 
the  mean  (which  is  0)  and  covariance  R  —  [eg^]  of  the  normally  distributed 
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range  measurement  noise  (e) .  v  -  WS(0,Pq)  defines  the  mean  (which  is  0)  and 
covariance  (Pg)  o£  v  in  the  weaker  form  called  wide  sense  (WS)  which  does  not 
make  any  distributional  assumptions.  Pq  is  known  as  the  error  covariance 
matrix  and  is  given  by 

Po  -  E[vyT]  -  E[(fi  -  (5.3) 

Since  the  initial  conditions  are  actually  random  variables,  the  system  state 
(g)  cannot  be  propagated  with  certainty.  Therefore,  the  state  itself  must  be 
considered  a  random  variable  and  the  result  is  a  stochastic  process  [29]. 

Following  the  Goldberger-Theil  method  given  in  [25],  the  idea  is  to  use 
5o  as  a  measurement 
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Now,  solving  for  ^  yields  the  Bayesian  solution  (Ip) .  As  derived  in  reference 
[30],  the  Bayesian  solution,  which  incorporates  the  a  priori  information  into 
the  linear  model,  is  given  by 

Ip  =  [P5I  +  H'^R-1H]-i  [PS^So  +  H^R-ly]  (5.5) 

The  linear  recursive  dynamic  model  is  obtained  by  placing  the  estimate  |j^  in 
the  role  of  the  a  priori  information 

Next,  the  goal  is  to  rewrite  equation  5.5  in  a  recursive  form.  The 
recursive  state  model  is  given  by 

fik+i  “  Hk+i  (5-6) 

where  ♦  is  the  state  transition  matrix  which  represents  how  a  dynamical  system 
naturally  evolves  from  one  state  to  the  next  in  the  absence  of  a  driving 
function.  (Note  that  the  system  would  be  in  a  constant  state  if  ♦  equals  the 
identity  matrix.)  Also,  represents  the  process  noise  with  associated 

system  covariance  matrix,  Q  -  E [ww^l . 

Following  a  procedure  given  by  reference  [26],  equation  5.5  can  be 
written  as 


h  -  (Pk/k-i  +  +  ElK-%)  (5.7) 

where  is  the  predicted  estimate  of  based  upon  observations  up  to 

k-1,  and  Pjc/ic-i  “  ‘^o''^(Sk/k- 1  "  ^k^  •  Equation  5.7  is  the  recursive  form  of  the 
Bayesian  estimator. 

Equation  5.7  can  then  be  transformed  into  the  Kalman  filter  measurement 
update  equation  by  using  the  Matrix  Inversion  Lemma  |The  Matrix  Inversion 
Lemma  is  used  to  obtain  the  following  relation:  [P'^  +  H’^R'^H]'^  «  P  - 
PH^(HPH'^  +  R)-'HP  [18]) 

^k  *■  Ik/k-l  *(k(Xk  -  ^k^k/k-l)  (5.8) 
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where  Yj^  are  the  new  measurements,  ^re  the  predicted  measurements 

based  on  the  last  k-1  estimates,  and  Kj^  is  called  the  Kalman  gain,  ^k/k-l 
represents  the  state  estimate  before  the  new  measurements  are  available,  while 
represents  the  state  estimate  after  the  new  measurements  are  incorporated. 
The  term  (Yj^  -  j )  forms  what  are  called  the  innovations;  it  represents 

the  part  of  Yj^  that  cannot  be  predicted  from  the  previous  measurement  data. 

The  innovations  are  the  same  as  the  residuals  of  the  OLS  estimator  (sj^)  in  the 
sense  that  the  Kalman  filter  minimizes  the  error  in  the  estimate  by  making  the 
innovations  white  (i.e.  -ero  mean  and  random). 

Equation  5.8  may  be  interpreted  as  follows:  the  new  "best"  state 
estimate  equals  the  previous  "best"  state  estimate  plus  the  new  information 
that  is  received  (the  innovations)  multiplied  by  a  gain  factor.  The  state 
estimates  are  the  best  in  the  sense  that  the  sum  of  the  squares  of  the  errors 
are  minimized.  This  makes  the  fxlter  "matched"  to  the  measurement  data  being 
processed  by  it.  When  the  Kalman  filter  model  is  matched  to  the  data,  the 
best  linear  unbiased  estimator  (BLUE)  is  obtained  [18]. 

Following  reference  [25],  the  gain  factor  is  deterained  from 

Kk  “  Pk/k-lH5[(HuPk/k-lHj  +  R)-'  (5.9) 

And  the  error  covariance  is  updated  by 

Pj,  =  (I  -  KkHk)Pk/k-l  (5.10) 

The  error  covariance  is  propagated  by  the  following  equation 

PkM/k  =  fPk*^  +  Q  (5.11) 

where  Q  represents  the  uncertainty  in  the  state  model.  Figure  5.1,  which 
emulates  a  figure  given  in  reference  [31],  shows  how  equations  5.6,  5.8,  5.9, 
5.10  and  5.11  interact  with  each  other  to  form  the  Kalman  filter. 

Note  that  ♦,  Q,  and  R  are  not  time  varying  in  the  above  equations  (no 
follows  them).  This  may  not  always  be  true.  For  example,  the  system 
dynamic  model  (4)  may  need  to  change  with  time,  or  the  system  covariance 
matrix  (Q)  may  be  updated  by  some  type  of  adaptive  scheme. 

The  major  point  of  this  section  is  that  the  Kalman  filter  is  the  same  as 
an  OLS  estimator  made  into  a  recursive  process  by  combining  the  incoming 
measurement  data  with  some  a  priori  information. 


5.1.3  The  Extended  Kalman  Filter  (EKF) 

Since  the  equations  that  relate  the  measurements  to  the  state  vector  are 
usually  nonlinear  (i.e.  the  H  matrix  is  nonlinear),  an  extended  Kalman  filter 
(EKF)  is  needed.  Therefore,  a  linearization  procedure  is  performed  when 
deriving  the  Kalman  filter  equations.  This  is  shown  below. 

First,  the  nonlinear  discrete-time  system  model  is  as  follows 

Xk  -  gk(fik)  +  4B  +  e  (5.12) 
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Enter  the  a  priori  estimate  3q 
and  its  error  covariance  Pq 


T 


Figure  5 . 1 


Kalman  filter  equations  (from  reference  (31)). 
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where  gj^  is  a  nonlinear  vector-valued  function  which  relates  the  system  state 
to  the  measurements  (see  reference  [32]).  Then,  the  linearized  model  is  given 
by 


^8k(^k/k-  1 ) 

a^k 


+  AB  +  e 


(5.13) 


The  gj^(^j^y I )  _term  signifies  that  the  linearization  is  taking  place  around 
the  estimate,  ^k/k-l*  Incorporate  equation  5.12  into  equation  5.7  to  obtain 
the  recursive  update  equation  for  the  extended  Kalman  filter 

k  “  ik/k-i  +  (Pk/k-i  +  hTr-1Hi,)-1(HTr-16Yi,)  (5.14) 

where  is  the  difference  between  the  actual  measurements  and  the  predicted 

measurements  (the  innovations)  as  shown  below 


k  — k  ”  8k^^k/k-l^ 

And  as  given  by  reference  [33] 

^8k(^k/k-l) 

H  =  - 

a^k 


(5.15) 


(5.16) 


As  will  be  seen  in  Section  5.6,  the  extended  Kalman  filter’s  update 
equation  as  given  by  5.14  is  very  useful  when  incorporating  ridge  regression 
theory  into  the  Kalman  filter. 


5.1.4  Selection  of  the  Pq  and  Q  Matrices  in  the  Kalman  Filter 

An  optimum  unbiased  estimator  arises  when  both  the  model  and  the 
estimator  match  the  process  which  generates  the  data.  Kalman  filter 
optimization  techniques  include  the  selection  of  Pq  and  Q  based  upon  past 
experience  or  by  adaptively  tuning  the  filter  until  its  innovations 
(residuals)  become  white  (i.e.  zero  mean  and  random).  However,  the  selection 
of  the  proper  Q  matrix  is  usually  not  a  very  easy  or  straightforward  task 
[31]  . 


For  example,  the  system  covariance  matrix  Q  is  often  set  artificially 
high  such  that  the  Kalman  filter  can  track  the  vehicle  when  it  encounters 
dynamics  such  as  turning- induced  accelerations.  Therefore,  the  Kalman  filter 
allows  more  noise  in  the  solution  during  periods  of  low  dynamics.  Although 
some  methods  exist  for  selecting  Q  adaptively  [34],  these  are  stochastic  in 
nature. 

Another  problem  arises  when  the  Kalman  filter  is  subjected  to  a  poor 
geometry  condition.  In  the  case  of  Inaccurate  Pq  and  Q  matrices  (a 
mismatched  filter),  the  filter  may  become  biased.  The  performance  of  a  biased 
Kalman  filter  is  not  readily  understood,  as  the  Kalman  filter  is  optimal  and 
defined  in  a  Gaussian  environment  only.  Therefore,  the  performance  of  a 
Kalman  filter  for  a  deterministic  maneuver  in  a  poor  geometry  condition  cannot 
be  predicted  from  the  regular  Kalman  filter  equations.  Furthermore,  near 
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colllnearity  effects  cannot  be  minimized  by  observing  the  innovations  in  the 
Kalman  filter’s  update  equation  (equation  5.8),  because  these  effects  only 
appear  in  the  estimator  This  is  shown  in  reference  [35]. 

In  applying  the  Kalman  filter  to  navigation,  the  MSE  is  the  appropriate 
error  criterion  to  be  minimized.  The  MSE  of  a  mismatched  Kalman  filter  is  not 
necessarily  the  smallest  obtainable.  Recall  that  the  MSE  is  the  sum  of  two 
components;  the  variance  term  and  bias  term  squared.  Since  the  Kalman  filter 
is  developed  from  the  OLS  estimator,  it  is  inherently  an  unbiased  estimator. 
This  restricts  the  MSE  equation  to  only  one  component  -  the  variance  term. 

The  next  section  presents  a  biased  estimator  derived  from  ridge 
regression  theory.  A  biased  estimator  is  purposely  not  matched  to  the  process 
that  generates  the  data  in  order  to  achieve  a  smaller  MSE. 


5.2  Theory  of  Ridge  Regression 

This  section  covers  a  biased  estimation  technique  called  ridge 
regression.  First,  the  ridge  estimator  is  analyzed  with  respect  to  the 
ordinary  least  squares  (OLS)  estimator  through  a  comparison  of  mean  square 
error  (MSE)  equations.  Second,  the  effects  of  range  measurement  noise  and 
bias  errors  on  the  ridge  estimator’s  solution  are  addressed.  Finally,  a 
linearized  recursive  ridge  processor  is  developed  which  explains  the  behavior 
of  the  extended  Kalman  filter  (EKF) . 


5.2.1  Historical  Perspective 

The  theory  of  ridge  regression  was  developed  by  Arthur  Hoerl  and  Robert 
Kennard  in  1970  [36]  based  on  some  earlier  work  by  Hoerl  in  1959  on  what  he 
termed  "ridge"  analysis  because  the  results  he  obtained  formed  what  looked 
like  ridges  in  the  output  data  [37].  The  purpose  for  the  development  of  ridge 
regression  was  to  combat  the  ill-effects  of  near  collinearity  when  it  arises 
in  the  predictor  matrix  of  a  linear  regression  model.  In  order  to  counter 
situations  with  near  collinearity,  ridge  regression  introduces  a  biasing 
parameter,  k.  Statisticians  have  had  much  debate  over  the  proper  selection  of 
this  biasing  parameter  [38,39,40]  (also  see  Oman  [41]  for  a  confidence  bound 
approach  to  selecting  k),  let  alone  the  fact  that  many  feel  that  a  biasing 
parameter  should  not  be  added  at  all  (see  Efron  [42]  who  points  out  the  heated 
controversy  among  statisticians  between  biased  and  unbiased  estimation).  The 
biasing  parameter  in  a  sense  "unbalances"  the  normal  linear  regression  model 
and  causes  a  bias  to  exist  in  the  solution.  This  is  in  direct  conflict  with 
an  unbiased  estimator  like  the  OLS  estimator  which  minimizes  the  sum  of  the 
squares  of  the  errors.  (This  was  examined  in  Section  4.3.)  Further,  detailed 
information  about  ridge  regression  may  be  found  in  references  [43,44,45]. 

In  1988,  Robert  Kelly  introduced  a  ridge  regression  signal  processor  for 
navigation  applications  to  combat  the  effects  of  nearly  collinear  range 
measurements  [46].  Inaccuracies  in  the  position  solution  become  highly 
inflated  when  range  measurements  are  nearly  collinear.  The  term  poor  geometry 
is  used  for  cases  when  range  measurements  are  nearly  collinear.  In  subsequent 
developments  over  the  past  two  years  [30,47,48,49],  it  has  been  determined 
that  a  proper  selection  of  the  biasing  parameter  can  be  made  based  on  the 
geometry  of  the  range  measurements.  Therefore,  ridge  regression  may  be  used 
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to  improve  position  estimates  obtained  from  a  navigation  system  in  the 
presence  of  poor  geometry. 


5.2.2  Biased  Estimation 

Position  fixing  using  any  type  of  navigation  system  suffers  a  loss  of 
accuracy  when  the  range  measurements  are  nearly  collinear,  i.e  poor  geometry. 
Under  such  conditions,  the  range  measurement  errors  are  propagated  into 
position  errors  which  are  highly  inflated  (see  Section  A. 3.1).  To  improve 
estimation  properties  when  several  parameters  are  to  be  estimated,  the  ridge 
estimator  deliberately  induces  biases  in  the  solution.  This  results  in  a 
solution  with  a  consistently  smaller  mean  square  error  (MSE)  than  the  ordinary 
least  squares  (OLS)  solution. 

Recall  that  the  MSE  consists  of  the  siim  of  two  components;  MSE  = 

Variance  +  Bias^.  The  ridge  estimator  takes  advantage  of  an  extra  degree  of 
freedom,  the  bias  term,  which  is  not  used  by  a  OLS  estimator.  In  effect, 
small  biases  induced  by  the  ridge  estimator  decrease  the  variance  term  such 
that  the  overall  MSE  is  smaller  than  the  MSE  obtained  from  an  unbiased 
estimator,  as  Illustrated  in  Figure  5.2. 

For  navigation,  the  most  useful  error  criterion  is  the  MSE.  It 
expresses  the  deviation  of  the  vehicle  with  respect  to  its  intended  path.  The 
ridge  estimator  has  the  property  that  the  MSE  in  the  presence  of  poor  geometry 
is  much  smaller  than  the  MSE  of  a  conventional  OLS  estimator  as  discussed 
above.  It  should  be  noted,  however,  that  the  ridge  estimator  cannot  remove  a 
steady-state  bias  in  the  estimator  caused  by  measurement  bias  errors.  This  is 
shown  in  Section  5. 2. 5.1. 


5.2.3  The  Ridge  Estimator 

Following  Kelly  [A8],  who  extends  the  ridge  regression  concept  as 
developed  by  Hoerl  and  Kennard  [A4]  for  navigation  applications,  the  linear 
model  for  a  system  with  an  unknown  n  x  1  measurement  bias  vector,  AB,  and  a 
measurement  noise  vector,  e,  is  given  by 

Y  -  M  +  e  (5.17) 

where  Y  is  the  n  x  1  range  measurement  vector,  §,  is  the  p  x  1  \inknown  system 
state  vector  (or  parameter  vector),  and  H  is  the  n  x  p  predictor  (or  design) 
matrix  which  relates  the  range  measurement  vector  to  the  system  state  vector. 
Also,  the  measurement  noise  is  uncorrelated;  cov(e  +  ^)  ••  [ee^]  ■  o^I,  where 
I  is  the  n  X  n  identity  matrix.  Recall  that  the  OLS  estimate  of  equation  5.17 
is 


loLS  “  (H'^H)-!  H'^Y  (5.18) 

The  corresponding  ridge  estimate  of  equation  5.17  is 

-  (H^'H  +  Pg)*^  H’^’Y  (5.19) 

where  Pj^  is  the  ridge  parameter  matrix  (which  is  formed  from  the  chosen  k 
parameters  as  will  be  seen  in  Section  5. 2. A).  When  P^  consists  of  zeros,  the 
ridge  estimator  reduces  to  the  OLS  estimator.  Adding  a  non-zero  ridge 


32 


ORX>C\^\RY  LEAST  SQUARES  (OLS) 


TRUTH ( p ) 


RIDGE  REGRESSION 

TRUTH ( p ) 


MSE  =  BLAS^  +  VARIANCE 


G?X  C29  /  0C6  EJ8 


Figure  5.2  Comparison  of  unbiased  and  biased  solutions  and 

definition  of  the  MSE. 
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parameter  matrix  to  purposely  upsets  the  balance  between  the  first  and 
second  moment  components  of  thereby  introducing  a  bias.  Substituting 
equation  5.17  into  equation  5.19  gives 

Ir  =  +  G^H^e  (5.20) 

where  G^  «  (H^H  +  Pr)'^.  From  equation  5.20  the  expected  value  (or  mean)  of 
the  ridge  estimator  is 

E[1r]  “  G^H^'Hg  +  G^H^'AB  (5.21) 

The  following  relation  is  introduced 

-  G^Pr^  (5.22) 

which  will  be  used  to  rewrite  equation  5.21.  Equation  5.22  can  be  proven  as 
follows 

(G^H^H)5  -  (I  -  G^Pr)^ 

-  G^'^G^Pr 
-  G,-l  -  Pr 
G^  -  (H'^H  +  Pr)*^ 

Therefore,  the  expected  value  given  by  equation  5.21  may  be  rewritten  as 

EL^r]  -  -  G^Pr^  +  G^H^AB  (5.23) 

Then,  the  ridge  estimator’s  solution  bias  is  seen  to  be 

BIASI^r]  -  -G^?r5  +  G^H'^4B  (5.24) 

The  variance  of  the  ridge  estimator  is 

VAR[iR)  -  OUTRAGE [G^H'^H]  (5.25) 

where  the  TRACE  operator  computes  the  sum  of  the  diagonal  terms  of  a  matrix. 

It  can  be  shown  that  the  bias  term  is  a  monotonically  increasing 
function  of  k  while  the  variance  term  is  a  monotonically  decreasing  function 
of  K  (see  reference  [35]).  This  implies  that  a  k  value  exists  which  will  give 
a  minimum  MSE  for  the  ridge  estimator.  (Recall  that  the  MSE  is  defined  as 
the  sum  of  the  variance  term  and  the  bias  term  squared.)  Similar  to  the 
derivation  of  the  MSE  for  the  OLS  estimator  in  Section  4.5,  the  MSE  of  ^r  is 
given  by 

MSE[fiR]  -  E[(1r  -  fi)T(lR  -  fi)] 

MSELIr]  -  E((-G^PRa  +  G^H^^  +  G^H’‘e)’'(-G^PRg  +  G^H^AB  +  G^H^'e)] 
mseiIr]  -  fi^pjGjG^PRg  +  m'^hgJg^h’^ab  -  ^’'pJgJg^h’'^  -  M’^HG^G^PrS 
+  OUTRAGE  [G^G^H’^'H] 

Since  G^  and  Pr  are  symmetric  matrices,  their  inner  products  are  squares  (i.e. 
gJg^  -  G^  and  PrPr  ■  Pr)*  This  yields  the  following  equation  for  the  MSE 

MSEliR]  -  -  25’'G^PrH'''M  +  M’’HG^H^AB  +  OUTRAGE [G^H'^H]  (5.26) 
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Reference  [35]  proves  that  the  MSE[4r)  is  less  than  or  equal  to  the  MSE[£ol5] 
based  on  a  "slope  argument":  The  slope  of  the  above  MSE  equation  can  be  shown 
to  be  negative  at  k  “  0.  This  implies  that  as  k  becomes  greater  than  zero, 
the  MSE  will  decrease. 


5.2.4  Selection  of  the  Ridge  Biasing  Parameter  Based  on  Geometry  Conditions 

The  problem  is  to  determine  the  range  of  values  for  which  the  MSE  of 
the  ridge  estimator  is  smaller  than  the  MSE  of  the  OLS  estimator  [50] 

MSE[4r]  -  MSE[1ols]  <0  (5-27) 

The  idea  is  to  transform  equation  5.27  into  its  canonical  coordinates 
such  that  the  proper  selection  of  the  parameters  can  clearly  be  seen. 
Canonical  coordinates  are  obtained  by  implementing  a  singular  value 
decomposition  (SVD)  of  H 

gnxp  .  unxn  j,nxp  (yT^pxp  (5.28) 

where  U  and  are  defined  as  orthogonal  transformation  matrices  that 
diagonalize  H  (i.e.  make  the  measurements  independent),  and  S  is  a  diagonal 
matrix  whose  elements  are  the  singular  values  of  H,  which  are  given  by 

V^.  Also,  where  is  a  p  x  p  diagonal  matrix  with  elements 

(the  eigenvalues).  Additionally,  the  relations  a  -  and  Aa  -  U^AB  exist 

where  V”  and  are  used  to  transform  the  system  state  vector  and  measurement 
bias  vector  into  their  corresponding  canonical  vectors.  The  ridge  parameter 
matrix  is  also  expressed  in  terms  of  its  canonical  coordinates  given  by 


Pr  -  vTpj^V 


Kj 

0 


(5.29) 


where  P^  is  a  p  x  p  diagonal  matrix  which  permits  a  different  ridge  biasing 
parameter  (kj^)  to  be  selected  for  each  element  of  the  canonical  system  state 
vector  a. 

In  order  to  perform  the  transformation  and  rewrite  equation  5.27  in 
terms  of  its  canonical  coordinates,  first  replace  B  with  S,  g  with  Va,  ^  with 
UAa,  and  Pj  with  P|  in  the  ridge  estimator’s  MSE  equation  (equation  5.26) 

MSE[gg]  -  (Vg)T(2'^S  +  P|)-^P^^(Va)  -  2(Vg)^(z’'z  +  P|) ' ^P|2£'r(UAg) 

+  (UAg)’'r(L^r  +  P|)‘2r’'(UAg)  +  o^TRACE[  (r'^s  +  P|)-22'^X] 

Since  2  is  diagonal,  2^2  -  2^  and  2^  “  2.  Since  V  and  U  are  orthogonal 
matrices,  then  V^V  -  I  and  U^U  ■  I.  This  reduces  the  previous  equation  to 

MSE[gR]  -  a’^(22  +  P|)-2pc2a  -  20^(2^  +  P|)-^P|2Ag 

+  ^’'‘2(2^  +  Pg)-22^  a^TRACE[(2^  +  P|)‘^2^] 
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Now,  rewrite  the  above  equation  from  its  matrix  form  to  its  subscript  form, 
given  the  fact  that  the  diagonal  elements  of  2^  are  given  by  and  P|  are 
given  by  for  i“l,p 


MSEI^r]  -  £ 
i-1 


(Xi  +  (Xj^  +  (X^  +  (Xj_  +  K^) 


(5.30) 


Similarly,  perform  the  same  manipulation  on  the  MSE  of  lo, e  given  in  Section 
A. 5 

MSE[ioLsl  “  Ab’'H(H'^H)-2h'^AB  +  OUTRAGE  [  (H^H)  *  M 
mselIols^  “  (U4a)‘^2(r’'z)'^2:’‘(U^)  +  outrage [ (s^'r) -^] 

MSE[ioLsl  ■  Aa’'z:(s2)-2rAa  +  OUTRAGE  [  (2^ )  ‘  M 
MSE[1ols]  -  Aa'^S'^Aa  +  OUTRAGE [2* 


MSE[^qls]  “  2 
i-1 


Aa^ 


(5.31) 


Now,  combining  equations  5.30  and  5.31  into  inequality  5.27  gives 


P 

2 

i-1 


0^X. 


Kj_ai 


+  ~- 


Aa|Xi 


(Xj  +  Kj^)^  (Xj  +  K^)^  (Xj^  +  Kj^)' 


2k  AaJ 

(Xj  +  Kj)^  Xj  Xj 


<  0 


(5.32) 


Inequality  5.32  will  be  satisfied  for  each  i  if 


o^Xj^  (Kj^a^-Aoj^V^)  ^  Aa| 


(Xj  +  Kj^)^  (Xj^  +  Kj)' 


(5.33) 


Simplifying  equation  5.33 
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<  (Xj^  +  Kj)^Aaj  +  (X^^  +  K^)^)a^ 

a^k\  +  Xj^(K^a|  -  2k j^Oj^VX^Aoj^  +  Aa^Xj)  -  (X^  +  2kj^Xj^  +  K|)Aa| 

-  (X^  +  2k^Xj^  +  <  0 

a^\\  +  t<|a|Xj^  -  2k j^aj^-A^Xj^Aaj^  +  Aa^X^  -  Aa|x|  -  Zkj^Xj^Ao^  -  K|Aa| 

-  o^Xj  -  Zkj^Xj^o^  -  k|o^  <  0 

K^a^X^  -  ZK^aj^-A^Xj^Aoj^  -  2K^X^Aa|  -  k^Ao^  -  Zkj^X^o^  -  k^o^  <  0 


Ki(Xia^  -  Aa^ 


a^)  -  ZCXj^a^AaiVXj  +  X^(Aa^  +  a^))  <  0 


Therefore, 


i 


3i I K ^  ^ 


Thus  for  K^  >  0 


1 


(5.34) 


1.  If  32  <  0,  aj  >  0,  there  is  no  Kj  >  0  which  satisfies  5.34 
and  Kj^  ■  0,  the  value  for  the  OLS  estimate. 

2.  If  aj  >  0,  aj  <  0,  then  5.34  is  satisfied  for  all  Kj^  >  0. 

3.  If  32  <  0,  aj  <  0,  then  a  "large"  Kj^  satisfies  5.34. 

4.  If  32  >  0,  aj  >  0,  then  a  "small"  Kj^  satisfies  5.34. 


These  restrictions  on  must  be  satisfied  if  the  ridge  estimator  is  to  have  a 
smaller  MSE  than  the  OLS  estimator  whenever  a  bias  component  AB,  as  well  as  a 
variance  component  o^  is  inflated  by  near  collinearity . 

In  practice,  the  vector  a  ■  V^g,  which  is  unknown,  is  replaced  by  its 
estimate  Sols  "  equation  5.34.  Also  required  is  some  information 

about  the  vector  Aa  -  U^AB.  Normally,  for  the  position  fixing  problem,  only 
nominal  values  of  Aoj  are  necessary  because  the  selection  of  Kj^  using  5.34  is 
not  a  strong  function  of  Aa^^  [49]  . 


A  procedure  has  been  developed  to  properly  choose  the  values  of  k^  based 
on  the  geometry  condition.  The  geometry  condition  is  represented  by  the 
eigenvalues  (X.)  of  H  H.  Poor  geometry  exists  when  the  eigenvalues  are  small. 
As  the  eigenvalues  become  small,  a  biasing  parameter  (k^)  is  added.  This 
limits  the  minimum  value  that  the  eigenvalues  can  obtain.  But  another 
restriction  exists  on  the  selection  of  the  parameters.  This  is  called  a 
convergence  criterion. 
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5.2.5  Selection  of  the  Ridge  Biasing  Parameter  Based  on  a  Convergence 
Criterion 

The  selection  of  the  ridge  biasing  parameter  also  depends  on  a 
convergence  criterion.  This  convergence  criterion  becomes  important  when  bias 
errors  exist  on  the  range  measurements  or  when  there  is  a  large  offset  in  the 
initial  estimate  of  position. 


5.2.5. 1  Convergence  Criterion  when  Range  Measurement  Bias  Errors  Exist 

Given  the  fact  that  a  poor  geometry  condition  exists,  the  ridge 
estimator  instantaneously  reduces  the  variance  inflation  when  noise  errors 
exist  on  the  range  measurements.  Additionally,  when  bias  errors  are  on  the 
range  measurements,  the  ridge  estimator  initially  reduces  the  bias  inflation, 
but  then  grows  with  a  certain  time  constant  of  convergence  to  the  OLS  inflated 
solution  bias.  This  can  be  seen  by  performing  the  following  analysis. 

Recall  that  the  bias  of  the  ridge  estimator  is  given  by 

BIASiIr]  -  +  pR)*^PRfi  +  (H’'H  +  Pr)-^h’‘AB  (5.35) 

Transforming  this  term  into  its  canonical  coordinates  yields 


BIAS[lR]i 


vq 


*^1 


aj  + 


*^1 


(5.36) 


The  second  term  of  equation  5.36  counteracts  the  inflation  of  the  measurement 
bias  errors  caused  by  poor  geometry.  By  inserting  a  positive  value,  there 
is  an  initial  shrinking  of  the  solution  bias.  This  can  be  seen  more  clearly 
by  looking  at  the  following  example. 


For  -  0.012; 


OLS  solution  (k^  “  0); 


-  AOi  -  9.13Aai 


vq 

Ridge  solution  with  -  0.05;  -  -  1.76^0^ 

Xi  + 

About  a  factor  of  5  improvement  in  the  solution  bias  can  be  expected  initially 
by  using  the  ridge  estimator  in  this  example.  However,  the  ridge  solution 
bias  exponentially  converges  to  the  OLS  solution  bias  as  shown  below. 

Start  with  equation  5.36.  Replace  the  two  coefficients  in  equation  5.36 
with  C  and  D 
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BlAS[gR]i 


(5.37) 


where 


-  -COi  + 


Xi  +  Kj_ 


VXi 

D  -  - 

Xi  +  Ki 


Now,  assume  only  measurement  bias  errors  exist,  and  propagate  the  bias  errors 
(ACi)  into  the  position  solution  (Oi).  Note  that  as  a  consequence  of  the 
ridge  estimator,  truth  ■  -  estimate  (see  equation  5.24  or  5.36). 

Perform  iterations  (n) ; 


n-1 : 

Oi  -  0, 

BIASCIrJi 

n-2; 

“  -DAa^, 

BIAS[iR]i 

n-3 ; 

-  -(C+l)DAaj^, 

BlAStij^li 

n-*®: 

biaslIr].  - 

D 

-  Aa. 

1  -  C 


-C(0)  +  DAtti  -  DAOi 
-C(-DAai)  +  DAOi  ■  (C+l)DAai 

-C(-(C+l)DAai)  +  DAOi  -  (C^+C+DDAOi 


(5.38) 


where  the  (1  -  C)'^  term  is  the  result  of  the  geometric  series  given  by  (C“‘^ 
+  +  •••  +  1)  as  n  goes  to  ®.  Now  substituting  the  canonical  terms  back 
in  for  C  and  D  in  equation  5.38  yields 


BIAS  [Sr]  i - - - Afli 

1 - 

Xi  +  Ki 


1 

BlAS[lj^]i  “  -  ^1  ■*  "final  value"  (5.39) 

VTi 

Note  that  Ki  drops  out.  Therefore,  when  measurement  bias  errors  exist,  the 
ridge  solution  bias  converges  to  the  OLS  solution  bias  as  the  n^uDber  of 
iterations  goes  to  ».  The  associated  time  constant  of  convergence  will  be 
derived  below. 

Again  looking  at  the  iterative  process  that  formed  equation  5.38,  it  can 
be  seen  that  the  ridge  estimator’s  solution  bias  exponentially  increases  to 
the  OLS  estimator’s  bias.  This  is  portrayed  in  Figure  5.3.  Figure  5.3  shows 
tkat  after  1  iteration  (In),  the  ridge  solution  bias  converges  to  within  about 
37Z  of  the  OLS  solution  bias.  (After  2  iterations  the  ridge  solution  bias  is 
within  14Z  of  the  OLS  solution  bias.)  For  the  ridge  estimator,  a  time 
constant  (r)  can  be  defined  as  the  time  it  takes  to  converge  within  about  e'^ 
or  37Z  of  the  "final  value",  which  in  this  case  ij  the  OLS  solution  bias.  To 
find  the  time  constant  of  convergence,  consider  iteration  n 
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Figure  5. 


BIAS 


(l-e')  =  0.63 
(l-e'^)  =  0.86 
(l-e^)  =  1.00 


crx  029  /  008  E38 


3  Exponential  increase  of  the  ridge  estimator’s  solution 
bias  to  the  OLS  estimator’s  solution  bias. 


B1AS[2r]i  - 


(5.40) 


1  -  C“ 

-  DAfli 

1  -  C 

Let  the  BIAS[gjj]j^  given  by  equation  5.40  be  called  an  "initial  value". 
Determine  for  which  n  will 

"initial  value" 

-  -  1  -  e'^  -  0.63  (5.41) 

"final  value" 

r  1 

I  I 

I  1  -  C“  I 

I  I 

I  -  DAa^  1 

1  I 

I  1  -  C  I 

I  1 

L  J 

-  -  0.63 

r  1 

I  1  I 

I  1 

I  -  DAfli  I 

I  1 

I  1  -  C  1 

I  I 

L  J 

1  -  C°  -  0.63 
C°  -  0.37 
lnIC“)  -  lnIO.37) 
n(ln[C])  -  -1 

n--(ln(C])-^  (5.42) 

Replacing  C  in  equation  5.42  with  the  canonical  term  it  represents 


n 


In 


\l  +  Kj_ 


(5.43) 


Then,  the  time  constant  of  the  ridge  solution  bias  convergence  to  the  OLS 
solution  is  given  by 
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^BC  “  ■^’^ra 


-In 


(5.44) 


where  AT^  is  the  update  interval  time.  As  a  special  note,  reference  [35] 
shows  that  even  increasing  the  number  of  range  measurements  available  for 
determining  a  position  fix  will  not  decrease  the  solution  bias  inflation 
caused  by  poor  geometry. 

5 .2. 5. 2  Convergence  Criterion  when  a  Large  Initial  Offset  Exists 

In  order  to  define  a  convergence  criterion  for  the  selection  of 
values  when  a  large  initial  offset  Is  present,  again  consider  the  bias  term  of 
the  ridge  estimator 

BIASI^r]  -  +  Pr)*^Pi{1  +  (H^H  +  Pjj)-1h'^AB  (5.45) 

To  simplify  the  analysis,  assxame  that  no  measurement  bias  errors  exist,  (AB  = 
0),  therefore 

B1AS[1r]  -  -(H^H  +  PR)'*PRfi  (5.46) 

Transforming  this  term  into  its  canonical  coordinates  yields 

BIAS[Ir]i  =  .  -  Oj  (5.47) 

Xl  + 

Based  on  equation  5.47,  a  system  response  time  constant  (^s}^}  may  be  defined 
by  implementing  a  similar  procedure  to  the  one  given  in  Section  5. 2. 5.1.  The 
result  is  Identical. 


5.2.6  Ridge  Recursive  Filter  Development 

So  far,  ridge  regression  has  been  presented  using  batch  processing, 
whose  exposition  is  easier  to  follow  than  a  recursive  presentation.  Section 
5.1.1  points  out  that  the  performance  of  the  batch  and  recursive  presentation 
is  equivalent.  A  procedure  for  using  ridge  regression  to  explain  the  behavior 
of  the  Kalman  filter  will  now  be  discussed. 

Similar  to  the  derivation  of  the  recursive  Kalman  filter  in  Section 
5.1.2,  obtain  the  state  update  equation  (equation  5.7)  in  ridge  form  [51] 

^Rk  -  (Pw'k-I  +  PRk  +  HTR-X>*'(Pk/k-iiRk/k-l  ^  hTr-%)  (5.48) 

where  the  ridge  parameter  matrix,  Ppj^,  has  been  added  to  counteract  the 
effects  of  an  ill-conditioned  (Hj[R'*Hj^)  term  (which  is  caused  by  poor 
geometry) . 

Extending  the  work  done  by  Kelly  [51],  equation  5.48  can  be  linearized 
using  equation  5.13  from  Section  5.1.3  to  obtain  the  linearized  state  update 
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equation 


'  ^Rk/k-l  +  (Pk/k-1  PRk  hJr->H,,)-1(HTr-16Yi,)  (5.A9) 

where  Pj^j^  has  been  added  to  the  (P^/k-i  ■*■  term  in  order  to  reduce  the 

effects  of  near  collinearity.  Equation  5.49  forms  the  basis  for  the  extended 
Ridge/Kalman  filter.  Also,  the  update  equation  for  the  error  covariance 
matrix  is  given  by 

Pk  =  (Pk/k-i  +  PRk  +  H^R-iHi,)-!  (5.50) 

Two  cases  may  now  be  defined  in  which  ridge  regression  can  explain  the 
behavior  of  the  Kalman  filter.  First,  in  the  absence  of  system  dynamics,  the 
ridge  parameter  matrix,  Prjj,  is  functionally  equivalent  to  the  state  error 
covariance  matrix,  Pk/k-i*  Therefore,  the  EKF  has  similar  convergence 
properties  (equation  5.44)  when  the  model  incorrectly  represents  the  system 
and  Pk/K-l  small.  Second,  if  dynamics  exist,  Pj^^  related  to  both  Pk/k-l 

and  the  system  error  covariance  matrix,  Q,  given  by  the  following  relation 

(^Pi,*^  +  Q)-l  -  (*Pi,*^  +  Q)-^  +  Pri,  (5.51) 

where  in  this  case,  the  -  symbol  means  "relating  to".  Again,  the  EKF  has 
similar  convergence  properties  when  the  model  incorrectly  represents  the 
system  and  (iPj^*^  +  Q)*^  is  small. 

Usually  Pj^  is  small,  which  means  that  the  Pk/k-l  equation  5.49 

is  large;  therefore,  there  is  no  collinearity  problem.  Note  that  Pk  varies  as 
the  Kalman  filter  is  updated,  but  normally  a  constant  Q  is  added  which 
"limits"  it  (i.e.  puts  uncertainty  in  the  model).  As  seen  in  the  above 
equation,  Pj^  will  be  large  when  Q  is  large.  Q  is  chosen  large  for  dynamic 
situations  (i.e.  turn-induced  accelerations).  When  a  poor  geometry  condition 
exists  in  addition  to  the  dynamic  situation,  Ppk  added  to  counteract 

the  large  Q.  Therefore,  a  proper  P^k  can  be  chosen  to  incorporate  ridge 
regression  into  the  Kalman  filter. 

The  procedure  recommended  for  implementing  ridge  regression  theory  into 
the  Kalman  filter  depends  upon  three  different  situations  [49]: 

(1)  The  GDOP  is  small,  therefore  H([R'*Hk  is  well-conditioned.  Qk/k-l 
chosen  to  match  the  model’s  process  noise  covariance  in  the  usual  way.  gy  is 
unbiased . 

(2)  The  GDOP  is  large,  therefore  Hj^R'^Hk  is  ill-conditioned.  But 

(Pk|k-1  +  HvR'^Hk)  is  well-conditioned.  Do  nothing.  The  estimator  is  matched 
to  the  model  and  gy  is  unbiased. 

(3)  The  GDOP  is  large  and  H^R'^Hk  is  ill-conditioned,  ^k/k-l  identified 
in  the  usual  way  and  (Pk/k.i  S^R'^Hk)  is  also  ill-conditioned  because  the 
elements  of  Pk/k-i  are  not  large  enough  to  remove  the  near  collinearity.  The 
ridge  parameter  matrix  »  ^Rk»  is  chosen  to  reduce  the  near  collinearity  of 
(Pk/k-i  ^kR'^Hk)  using  the  selection  rules  given  by  equation  5.34. 

Situations  (1)  and  (2)  are  optimized  in  the  usual  way,  i.e.  choosing 
Pk/k-i  in  the  filter  such  that  the  innovations  are  white.  Situation  (3),  on 
the  other  hand,  minimizes  the  MSE  by  applying  equation  5.34  to  determine  P^k* 
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When  situation  (3)  occurs,  only  the  values  of  are  adjusted  by  adding 

P{5j^.  In  summary,  the  optimization  procedure  for  situation  (3)  is: 

1.  Identify  and  determine  Iqls* 

2.  Calculate  ^oLS  “  ^  “  U^AB.  Note  that  since  the 

measurement  bias  error  AB  is  unknown,  it  must  be  approximated. 

3.  Insert  the  calculation  results  from  step  2  into  equation  5.34 
to  obtain  for  i  »  1,2,  ...  p.  (This  gives  P|j^.) 

4.  Using  V^,  form  Pj^v  from  P^j^  using  equation  5.34  and  add  Pj^j^ 

to  (P^/k-l  make  it  well-conditioned. 

The  key  idea  in  developing  the  ridge  recursive  filter  is  the  following 
Each  step  in  the  recursive  process  is  viewed  as  a  new  prior  linear  model 
wherein  the  last  estimate  ^Rk/k-l  prior  equation  for  the  next 

iteration.  The  ridge  solution  is  recomputed  at  each  step  using  the  above 
selection  rules  to  determine  a  proper  Prr. 
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6. 


MULTISENSOR  NAVIGATION  SOLUTION  AND  INTEGRITY 


6.1  Multisensor  Navigation  Solution 

Navigation  measurement  data  can  be  expressed  in  terms  of  range  or 
bearing  with  respect  to  a  reference  station,  as  discussed  in  Chapter  2.  The 
relation  between  the  measurement  data  and  the  user  position  is  nonlinear,  as 
shown  by  the  following  two  meisurement  equations: 

Range:  R^  -  ((X-Xi)^  +  (6.1) 

Bearing:  6^  “  arctan{ (X-Xj^)  /  (Y-Yj^) }  (6.2) 

where  (X,Y)  is  the  user  position,  (Xj^,Yj)  is  the  location  of  reference  station 
i,  Rj^  is  the  range  between  the  user  and  reference  station  i,  and  ©.  is  the 
bearing  of  the  user  with  respect  to  reference  station  i  measured  clockwise 
with  respect  to  North,  see  also  Figure  2.2.  Another  widely  used  measurement 
is  the  time  difference  (TD)  between  the  time  of  arrivals  of  signals  from  two 
different  reference  stations.  The  TD  defines  a  hyperbolic  line-of-position 
and  the  measurement  equation  is  given  by: 

TD:  TDij  -  b  -  Rj  +  Rj  (6.3) 

where  b  is  the  length  of  the  baseline  between  station  i  and  j,  and  R^^,  R.  are 

given  by  equation  6.1. 

To  obtain  a  position  solutions,  equations  6.1  through  6.3  are  usually 
linearized  with  respect  to  some  reference  point  (X,^),  the  a  priori  estimate. 
Following  the  linearization  procedure  presented  in  Section  4.4.1,  equations 
6.1  through  6.3  can  be  linearized  as  follows. 


(6.4) 


(6.5) 


(6.6) 


Equations  6.4  through  6.6  relate  a  change  in  the  user  position  to  a  ch;’','’,  ^  in 
the  range,  bearing,  and  time  difference  measurements,  respectively.  In 
general,  each  of  these  equations  is  of  the  form: 

6j^  -  hi  fdX  6YI’’  (6.7) 

where  h^  is  the  row  vector  corresponding  to  measurement  number  i,  given  by  y^. 
Equation  6.7  can  be  written  to  include  all  thr  different  measurements  as 
follows . 


6Y  -  E  6Q, 


(6.8) 
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where  Y  is  the  measurement  vector,  Q,  is  the  user  state  vector,  and  H  is  the 
matrix  containing  the  row  vectors  given  by  h^.  Equation  6.8  can  be  used 
directly  in  a  recursive  estimator  such  as  the  least  squares  estimator  or  an 
extended  Kalman  filter,  as  discussed  in  Section  5.1.3. 

The  presence  of  an  unknown  variable  in  either  of  the  above  measurement 
equations  adds  another  variable  to  the  user  state  vector  For  instance,  if 
a  range  measurement  contains  an  unknown  clock  phase  offset,  then  the 
corresponding  pseudorange  measurement  equation  is  given  by: 

Pseudorange:  -  ((X-Xi)^  +  +  B  (6.9) 


where  B  is  the  unknown  clock  phase  offset.  The  linearized  measurement 
equation  is  then  given  by: 


6PRi 


(6.10) 


Similar  expansions  can  be  developed  for  other  unknowns,  e.g.  velocity  and 
acceleration,  which  are  to  be  determined  simultaneously  with  the  position 
solution.  This  concludes  the  development  of  a  unified  solution  for  a 
multisensor  position  solution. 


6.2  Multisensor  Integrity 

A  multisensor  navigation  system  has  many  redundant  measurements 
available  to  perform  fault  detection  and  isolation  (FDI),  which  is  also 
referred  to  as  receiver  autonomous  integrity  monitoring  (RAIM) .  The  basic 
idea  is  to  use  the  inconsistency  in  the  measurement  data  to  derive  failure 
detection  and  isolation  parameters.  The  next  section  presents  a  derivation  of 
the  least  squares  residual  vector  which  is  commonly  used  for  FDI.  The  least 
squares  residual  method  does  not  rely  on  the  history  of  the  measurement  data; 
it  is  based  on  a  "snapshot"  least  squares  solution.  The  use  of  this  method  is 
justified  for  several  reasons: 

1.  The  main  purpose  of  the  least  squares  residual  method  is  to  detect 
slowly  growing  measurement  errors  which  do  not  grow  rapidly  enough  to  be 
detected  by  an  input  data  editor.  (An  input  data  editor  would  simply 
reject  measurement  data  which  gives  rise  to  large  residuals). 

2.  Filtering  of  the  measurement  data  could  reduce  the  false  alarm  rate; 
however,  for  GPS,  the  measurement  errors  are  dominated  by  Selective 
Availability  which  has  a  correlation  time  on  the  order  of  minutes. 

Also,  the  least  squares  solution  still  allows  for  some  variance 
reduction  by  averaging  over  several  measurements. 
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3.  No  matter  what  kind  of  estimator  is  used,  the  position  bias  error 
caused  by  measurement  bias  errors  will  always  converge  to  the  position 
bias  error  of  the  least  squares  estimator  as  shovni  in  Chapter  5.  In 
almost  all  cases,  the  time  constant  of  convergence  is  short  compared  to 
the  slow  error  growth  of  a  difficult  to  detect  measurement  error. 

4.  Integrity  information  must  be  available  soon  after  receiver  start-up 
or  re-start. 

5.  The  integrity  performance  of  filtered  data  is  difficult  to  guarantee 
under  dynamic  conditions . 


6.2.1  The  Least  Squares  Residual  Vector 

With  reference  to  Section  4.2  and  Figure  4.1,  the  Projection  Theorem 
yields  the  ordinary  least  squares  (OLS)  estimator  given  by 

k  -  (HTH)-IhTy  -  ioLS  (6.11) 

Multiplying  equation  6.11  by  H  yields 

Hi  -  (6.12) 

X  -  PY  (6.13) 

where  the  projection  matrix  P  =  projects  the  measurement  vector  Y 

onto  the  estimation  space  defined  by  the  column  space  of  H.  Therefore,  the 
rank  of  P  for  an  overdetermined  system  is  equal  to  the  number  of  unknowns  in 
the  user  state  vector  jS.  This  projection  finds  the  closest  point  on  the 
estimation  space  with  respect  to  Y  by  constructing  a  perpendicular  line  from  Y 
to  the  column  space  of  H,  see  Figure  4.1. 

The  least  squares  residual  vector  ej^,  also  referred  to  as  the  error 
vector,  is  given  by 

eg  “  Y  -  Y  -  Y  -  PY  -  (I  -  P)Y  (6.14) 

where  (I  -  P)  is  also  a  projection  matrix,  it  projects  any  vector  Y  onto  the 
orthogonal  complement  of  the  estimation  space.  The  rank  of  (I  -  P)  is  equal 
to  the  degrees  of  freedom  of  the  overdetermined  system.  For  instance, 
consider  the  GPS  solution  which  requires  4  measurements.  If  6  measurements 
are  available,  then  the  rank  of  (I  -  P)  is  equal  to  2.  This  also  means  that 
the  number  of  rows  of  the  (I  -  P)  matrix  can  be  reduced  to  2. 

The  expected  value  of  e^  is 

£[6^]  -  (I  -  P)  E[Y]  -  (I  -  P)M  (6.15) 

where  AB  is  a  vector  containing  bias  errors  present  in  the  measurements.  The 
covariance  matrix  of  ej^  is 

COVCsk]  -  EteReJ]  -  E[(I  -  Dee^d  -  P)'^] 

-  (I  -  P)  COV[e]  (I  -  P)  (6.16) 

where  e  is  a  vector  which  represents  zero-mean  measurement  noise.  Note  that 


47 


(I  -  P)^  -  (I  -  P),  which  is  one  of  the  basic  properties  of  a  projection 
matrix.  The  second  property  of  a  projection  matrix  is  that  (I  -  P)^  *  (I  - 
P);  a  projection  followed  by  the  same  projection  does  not  change  the  result 
from  the  first  projection.  If  the  measurement  noise  is  uncorrelated  and 
normally  distributed  with  equal  variances,  then  the  covariance  of  the 
measurement  noise  is 

COV[e]  -  0^1  (6.17) 

where  a  is  the  standard  deviation  of  the  measurement  noise.  It  then  follows 
that  the  covariance  of  the  error  vector  is 

COV[e^)  -  0^(1  -  P)  (6.18) 

Using  the  Projection  Theorem,  the  error  vector  ej^  is  easily  obtained,  as 
given  by  equation  6.14.  This  error  vector  is  widely  used  for  failure 
detection  and  isolation,  see  for  instance  references  [52,  53,  54], 

Following  reference  [55],  the  above  general  solution  can  be  expanded  to 
include  a  weighting  matrix  W  to  represent  the  relative  importance  of  the 
measurements.  The  weighted  least  squares  estimate  is  obtained  by 

fiy  “  (H^WH)-^H’'WY  (6.19) 

where  W  is  a  positive  definite  weighting  matrix.  The  best  unbiased  weighting 
matrix  is  given  by  the  inverse  of  the  covariance  matrix  of  the  measurement 
vector  [55] 

W  =  (COV[Y] )*^  (6.20) 

Again,  the  projection  Hgy  onto  the  estimation  space  and  the  error  (Y  -  H^y) 
are  perpendicular.  It  then  follows  that  the  new  projection  matrix  is  given  by 

Py  =  H(H’'wH)'^H'’^W  (6.21) 

and  the  weighted  error  vector  is  given  by 

^RW  -  (I  -  Pw)I  (6.22) 

6.2.2  Fault  Detection  and  Isolation  Using  the  Least  Squares  Residual  Vector 

Once  the  least  squares  residual  vector  or  error  vector  has  been 
obtained,  several  methods  are  available  to  perform  the  failure  detection  and 
isolation.  References  by  Parkinson  and  Axelrad  [52]  and  by  Sturza  [53]  use 
the  square  of  the  magnitude  of  the  error  vector  as  the  basis  for  the  decision 
variable 

D  -  ejen  (6.23) 

If  the  measurements  are  normally  distributed  with  the  same  variance  o^,  then 
the  normalized  decision  variable 

-  D/o2  (6.24) 

has  a  chi-square  distribution  with  rank(I  -  P)  degrees  of  freedom.  If  the 
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measurements  have  nonzero  means,  then  has  a  noncentral  chi-square 
distribution.  The  noncentral  chi-square  distribution  and  calculation  methods 
for  this  distribution  are  documented  in  references  [52,  53]  . 

If  the  degrees  of  freedom,  m,  is  greater  than  one,  the  decision 
variable,  D,  could  be  scaled  by  the  degrees  of  freedom.  Reference  [52] 
suggests  the  use  of  the  range  residual  parameter  r  for  the  test  statistic  or 
integrity  parameter: 

r  »  V(D/m)  -  V(t|ej{/m)  (6.25) 

This  particular  integrity  parameter  has  also  been  adopted  for  the  baseline 
scheme  contained  in  the  document  "Minimum  Operational  Performance  Standards 
for  Supplemental  Airborne  Navigation  Equipment  Using  GPS"  [15].  Other 
integrity  parameters  have  been  suggested,  see  for  instance  references  [54,  56, 
57,  58,  59,  60]. 

Next,  the  integrity  parameter  must  be  compared  with  a  threshold,  T.  If 
the  threshold  is  exceeded,  a  fault  is  detected,  otherwise,  no  fault  is 
detected.  The  performance  of  the  integrity  algorithm  can  be  expressed  in 
terms  of  two  probabilities: 

Probability  of  a  false  alarm:  Pp^  “  P(  r>T  I  no  fault) 

Probability  of  a  missed  detection:  Pj^P  ••  P(  r<T  ]  fault) 

Chapter  3  lists  the  preliminary  requirements  for  supplemental  navigation  using 
GPS  for  these  probabilities:  Pr^  <  0.002/Hr;  and  P^^p  <  0,001  on  a  per  sample 
basis.  For  GPS,  it  is  generally  assumed  that  independent  samples  are 
available  every  2  minutes  (correlation  time  of  Selective  Availability).  It 
then  follows  that  the  Ppj^  <  0.000067  on  a  per  sample  basis.  From  reference 
[53],  the  threshold  T  must  be  set  at  approximately  V(16 .6a^)  for  the  case  of 
one  degree  of  freedom.  Assuming  a  standard  deviation  of  approximately  30 
meters,  the  threshold  T  ••  122  meters.  If  the  standard  deviation  of  the 
measurement  noise  is  100  meters,  the  threshold  T  “  407  meters. 

To  satisfy  the  probability  of  a  missed  detection  for  GPS,  it  follows 
from  reference  [53]  that  the  radial  protection  error  has  to  be  increased  to 
approximately  1500  meters,  rather  than  the  550  meters  required  for 
nonprecision  approaches.  It  should  be  noted  that  the  performance  of  the 
integrity  algorithm  would  improve  if  more  measurements  are  available. 

However,  not  all  these  measurements  are  required  for  positioning  and  FDI,  the 
question  remains  if  all  measurements  have  to  be  Included  in  the  FBI  process. 
Also,  other  detection  techniques  have  been  proposed  which  claim  to  perform 
better  than  the  above  method  [54].  Further  studies  are  recommended  to  analyze 
the  performance  of  several  integrity  techniques,  and  to  quantify  the 
performance  of  specific  multisensor  systems,  such  as  GPS/LORAN. 

In  addition  to  fault  detection,  several  techniques  have  been  proposed  to 
perform  fault  isolation,  such  as  the  maximum  likelihood  estimation  approach 

[53] ,  comparison  of  detection  thresholds  based  on  different  subsets  [52],  and 
coordinate  changes  to  maximize  the  visibility  of  a  certain  measurement  error 

[54] .  A  detailed  study  of  these  techniques  is  outside  the  scope  of  this 
report,  but  is  required  to  address  the  issue  of  sole  means  navigation. 
Therefore,  further  studies  are  also  recommended  in  this  area. 
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7. 


PROTOTYPE  HYBRID  GPS/LORAN  RECEIVER 


7.1  Hardware  Configuration 

The  block  diagram  of  the  hardware  configuration  for  the  prototype  hybrid 
GPS/LORAN  receiver  is  shown  in  Figure  7.1.  A  four-channel  GPS  receiver 
(Motorola,  model  Eagle)  and  an  eight-channel  LORAN-C  receiver  (Advanced 
Navigation,  Inc.,  Model  5300),  both  employing  continuous  tracking,  are  used  to 
collect  GPS  and  LORAN-C  data.  Only  the  raw  measurement  data  from  both 
receivers  is  used  to  determine  the  position  solution.  The  two  receivers  are 
interfaced  to  a  microcomputer  (model  AT)  through  two  serial  communication 
ports.  The  microcomputer  is  also  interfaced  to  a  course  deviation  indicator 
(CDI,  model  KI  206),  through  a  parallel  port,  to  display  the  guidance  data  to 
the  pilot.  All  of  the  hardware  used  is  commercially  available  equipment, 
except  for  the  interface  between  the  microcomputer  parallel  port  and  the  CDI 
instrument,  which  was  designed  and  implemented  at  Ohio  University. 


7.2  Software  Algorithms 

The  software  modules  implemented  on  the  hybrid  GPS/LORAN  receiver  are 
executed  in  realtime.  The  algorithm  is  given  by: 

initialization 
WHILE  in  operation 

DO  once  per  second 

check  for  keyboard  input  data 
IF  keyboard  input  data 

process  keyboard  data 

END 

check  for  GPS  and  LORAN-C  measurement  data,  and 
request  LORAN-C  data 
IF  sufficient  data 

calculate  position 
determine  integrity 

END 

update  CDI  and  status  screen 
store  all  relevant  data 

END 

END 

system  shut-down 


7.2.1  Initialization 

During  initialization,  the  following  tasks  are  performed: 

1.  Open  the  GPS/LORAN  initialization  data  file. 

2.  Read  the  following  data  from  the  initialization  data  file: 

a.  GPS  ephemerides 

b.  LORAN-C  transmitter  locations 

c.  LORAN-C  propagation  data 

d.  Waypoints 

e.  Position  estimate 

f.  Integrity  threshold 
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Figure  7.1  Hardware  configuration  of  the  prototype 
hybrid  GPS/LORAH  receiver. 
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g.  CDI  centering  and  offset  corrections 

h.  Ephemeris  threshold  (for  comparison  of  old  and 
new  ephemeris  data) 

3.  Close  the  GPS/LORAN  initialization  data  file. 

4.  Open  the  GPS/LORAN  output  file. 

5.  Install  serial  communication  interrupt  service  routines. 

6.  Output  the  status  screen. 


7.2.2  Check  For  and  Process  Keyboard  Input  Data 

The  keyboard  input  buffer  is  checked  for  input  data.  If  input  data  is 
present,  the  data  is  read  and  processed  according  to  the  list  of  possible 
commands  given  in  Table  7.1. 


Command 


Action 


1 

1  Q 

— 1 — 
1 

Initiate  receiver  shut-down 

- 1 

1 

1  Wxy 

1 

Waypoint,  from 

X  to  y,  where  x,y 

“0,1,... ,9.  1 

1  Fmx 

1 

Inject  failure 

into  signal  number  m;  x  is  the  failure  type  1 

( 

X  =  A,  a 

100  meter 

step 

failure  | 

1 

X  ■  B,  a 

1,000  meter 

step 

failure  | 

1 

X  ■  C,  a 

10,000  meter 

step 

failure  j 

1 

X  ■  D,  a 

1  m/s 

ramp 

failure  j 

1 

X  «  E,  a 

2  m/s 

ramp 

failure  | 

1 

X  =  F,  a 

3  m/s 

ramp 

failure  j 

1 

X  «  G,  a 

4  m/s 

ramp 

failure  j 

1 

X  «  H,  a 

5  m/s 

ramp 

failure  j 

1 

X  ■  I,  a 

10  m/s 

ramp 

failure  i 

1 

X  ■  J,  a 

25  m/s 

ramp 

failure  { 

1 

X  ■  K,  a 

50  m/s 

ramp 

failure  | 

1 

X  ■  L,  a 

100  m/s 

ramp 

failure  | 

1  R 

1 

Reset  injected 

failure  to  zero. 

1  A 

1 

Set  CDI  scale  to  ±1.25  nmi 

(approach  mode).  | 

1  E 

1  c 

1 _ 

1 

1 

_ 1 _ 

Set  CDI  scale  to  ±5.00  nmi 
Clear  the  display. 

(enroute  mode).  | 

1 

1 

Table  7.1  List  of  keyboard  commands  for  the  GPS/LORAN  receiver. 


7.2.3  Check  for  GPS  and  LORAN-C  Measurement  Data,  and  Request  LORAN-C  Data 

During  system  initialization,  the  GPS  receiver  is  commanded  to  send 
measurement  date  at  a  rate  of  once  per  second.  As  soon  as  GPS  data  is 
received,  a  LORAN-C  measurement  trigger  command  must  be  sent  to  ensure  that 
the  LORAN-C  data  is  valid  at  the  same  time  as  the  GPS  data.  Following  the 
measurement  trigger  command,  LORAN-C  data  is  requested  and  collected  for  up  to 
eight  receiver  channels.  All  data  is  verified  for  validity  as  indicated  by 
the  receivers.  If  at  least  five  measurements  are  valid,  sufficient  data  is 
available  for  the  position  calculation.  The  five  measurements  are  used  to 
solve  for  three-dimensional  position,  clock  offset  with  respect  to  GPS  time, 
and  clock  offset  with  respect  to  LORAN-C  time.  The  number  of  required 
measurements  could  be  reduced  to  four  if  the  GPS  and  LORAN-C  receivers  measure 
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the  time-of-arrivals  with  respect  to  the  same  clock,  and  if  the  hardware 
delays  of  both  receivers  are  known.  A  minimum  of  six  valid  measurements  are 
required  to  execute  the  integrity  algorithm. 


7.2.4  Calculate  Position 

The  hybrid  position  solution  is  based  on  a  least-squares  solution.  The 
GPS  and  LORAN-C  measurements  are  equally  weighted.  Reference  [4]  contains  a 
detailed  description  of  the  algorithm  used.  Since  the  measurements  from  GPS 
and  LORAN-C  are  equally  weighted,  the  accuracy  of  the  hybrid  system  will  be 
mostly  determined  by  the  LORAN  measurements.  For  this  effort,  standard  LORAN 
propagation  models  are  used  such  that  the  achieved  accuracies  are 
representative  for  current  LORAN  receivers.  Because  of  this,  the  accuracy  of 
the  hybrid  system  will  not  be  as  good  as  that  provided  by  GPS;  however,  the 
availability  and  integrity  of  the  hybrid  system  exceeds  that  of  GPS  by  several 
orders  of  magnitude  [4] .  At  the  same  time,  the  hybrid  system  accuracies  are 
still  well  within  the  current  requirements.  The  accuracy  of  the  hybrid 
GPS /LORAN  system  can  be  improved  upon  significantly  through  one  or  more  of  the 
following  methods: 

1.  Use  of  a  weighting  matrix  W  to  incorporate  the  statistical  knowledge 
of  the  measurements  (see  Section  6.1); 

2.  Calibration  of  LORAN  using  validated  GPS  positions; 

3.  Use  of  improved  LORAN  propagation  models  which  could  contain  seasonal 
correction  data  based  on  the  LORAN-C  monitor  network. 

Especially  the  latter  two  methods  are  very  promising,  these  methods  have  the 
potential  to  achieve  LORAN  measurement  accuracies  very  close  to  those  provided 
by  GPS  [4,  8,  61,  62,  63).  The  effects  on  the  integrity  performance  of  the 
first  method  is  not  well  understood  at  this  time.  Further  efforts  are 
recommended  with  respect  to  each  of  the  above  methods  to  improve  the  accuracy 
of  hybrid  GPS /LORAN. 


7.2.5  Determine  Integrity 

Integrity  is  calculated  using  the  range  residual  parameter  given  by 
equation  6.25  in  Chapter  6.  It  should  be  noted  that  the  performance  of  this 
particular  integrity  scheme  does  not  fully  demonstrate  the  integrity 
capabilities  of  a  hybrid  GPS/LORAN  system.  For  instance,  the  LORAN-C  system 
will  have  its  own  integrity  monitor  (also  known  as  aviation  blink),  which  is 
anticipated  to  be  fully  operational  by  1992.  At  that  time,  the  main  function 
of  the  receiver  autonomous  integrity  monitor  for  LORAN  would  be  to  detect  rare 
occurrences  of  receiver  cycle  slip.  One  of  the  purposes  of  this  effort, 
however,  is  to  evaluate  the  range  residual  parameter  technique,  which  is  the 
baseline  RAIM  scheme  proposed  in  reference  [15). 


7.2.6  Update  CDI  and  Status  Screen 

The  current  "from"  and  "to"  waypoints  are  used  to  calculate  the  course 
deviation  to  be  displayed  to  the  pilot.  The  CDI  needle  deviation  is  scaled 
according  to  the  phase  of  flight;  ±1.25  nmi  for  the  nonprecision  approach,  and 
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±5  nmi  for  enroute  navigation.  The  results  of  the  integrity  calculation  are 
used  to  drive  the  CDI  flag.  The  status  screen  is  used  to  monitor  the 
performance  of  the  GPS/LORAN  receiver,  and  to  observe  the  numerical  effects  of 
the  injection  of  signal  failures.  Parameters  of  interest  include;  current 
time,  position  estimate,  GPS  satellites  and  LORAN  transmitters  being  used, 
current  waypoints,  injected  failure  type  and  magnitude,  CDI  value,  integrity 
parameter  and  integrity  flag. 


7.2.7  Store  all  Relevant  Data 

All  data  collected  from  the  GPS  and  the  LORAN-C  receiver  is  recorded  on 
a  flexible  disk.  This  data  allows  for  a  complete  flight  evaluation  in  the 
laboratory  environment,  see  also  Chapter  8. 


7.2.8  System  Shut-Down 

During  system  shut-down,  the  following  tasks  are  performed; 

1.  Close  the  GPS/LORAN  output  file. 

2.  Save  all  new  ephemeris  data  collected  during  the  run. 

3.  Remove  serial  communication  interrupt  service  routines. 

7.3  Prototype  GPS/LORAN  Receiver  Installation 

The  prototype  GPS/LORAN  receiver  is  installed  in  a  Piper  Saratoga  PA- 
32-301,  N8238C,  which  is  owned  by  Ohio  University.  The  N8238C  is  a  1980  model 
aircraft  with  a  fixed  landing-gear,  and  a  useful  load  capacity  of  1,537 
pounds.  The  aircraft  is  equipped  as  a  flying  laboratory. 

The  GPS  antenna /pre-amplifier  is  mounted  on  top  of  the  fuselage  at  a 
distance  of  approximately  A  feet  from  the  front  windshield.  The  antenna/pre¬ 
amplifier  is  A. 5  inches  square,  and  2.5  inches  high.  Located  in  the  top  of 
the  pre-amplifier  is  a  micro-strip  antenna.  A  one-foot  slanted  LORAN  antenna 
is  also  mounted  on  top  of  the  fuselage,  approximately  8  feet  back  from  the  GPS 
antenna.  Both  antennas  are  connected  to  the  corresponding  receivers  which  are 
located  in  an  equipment  rack  together  with  the  microcomputer.  This  equipment 
rack  replaces  one  of  the  passenger  seats  in  the  back  of  the  airplane.  The 
microcomputer  is  also  connected  to  a  KI  206  course  deviation  indicator  mounted 
in  the  primary  view  of  the  pilot.  The  equipment  is  connected  to  the  aircraft 
power  system,  which  provides  lA  Volts  and  28  Volts  DC;  and  110  Volts  AC 
through  a  solid  state  power  Invertor.  The  prototype  receiver  can  be  installed 
in  less  than  one  hour. 
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8. 


POST-FLIGHT  REAL  TIME  GPS/LORAN  SIMULATOR 


A  post-flight  real  time  simulator  has  been  developed  to  support  the 
design  and  evaluation  of  the  GPS/LORAN  receiver.  The  simulator  consists  of  a 
microcomputer  which  simulates  the  GPS  and  LORAN  receivers  based  on  previously 
collected  flight  data.  Figure  8.1  shows  the  block  diagram  of  the  simulator 
configuration.  The  flight  data  is  read  by  the  simulator  computer  program  and 
transmitted  over  two  serial  ports  following  the  same  protocols  as  those 
employed  by  the  GPS  and  LORAN  receivers.  Initially,  the  simulator  used  data 
which  was  collected  in  September  of  1988  [4].  This  data  includes  several 
turns  and  data  outages.  Based  on  this  data,  the  real  time  software  could  be 
fully  exercised  in  the  laboratory  environment.  This  proved  to  be  a  very 
powerful  approach,  since  several  potential  problems  were  corrected  before  the 
actual  flight  tests.  These  problems  ranged  from  incorrect  ephemeris 
information  obtained  from  the  GPS  receiver  to  numerical  problems  due  to  the 
injection  of  large  signal  failures. 
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Fipure  8.1  Post-flight  real  time  GPS/LORAN  simulator  block  diagram. 


9. 


LABORATORY  TESTING  OF  THE  GPS/LORAN  RECEIVER  USING  FLIGHT  DATA 


GPS/LORAN  flight  data  collected  on  September  16,  1988  was  used  to 
evaluate  the  realtime  hybrid  GPS/LORAN  software.  This  particular  set  of 
flight  data  was  chosen  because  a  differential  GPS  (DGPS)  reference  trajectory 
is  available.  The  DGPS  trajectory  is  based  on  GPS  flight  data  corrected  for 
known  GPS  errors  as  determined  by  a  GPS  reference  station.  This  reference 
station  consists  of  a  GPS  receiver  and  processing  equipment.  The  accuracy  of 
the  DGPS  system  is  better  than  10  meters  (2  drms)  [4].  This  level  of  accuracy 
qualifies  DGPS  very  well  for  a  truth  reference  system  for  the  evaluation  of 
navigation  results  where  the  highest  accuracy  requirement  is  100  m  (2  drms). 
Reference  [4]  describes  the  DGPS  system  used  for  this  study  in  detail. 

Figure  9.1  shows  the  two-dimensional  position  errors  as  a  function  of 
time,  and  Figure  9.2  shows  the  corresponding  scatter  plot.  The  flight  lasted 
for  approximately  70  minutes.  Two  relatively  large  discontinuities  in  the 
data  are  caused  by  the  exchange  of  flexible  disks  and  by  a  system  restart. 

The  system  was  restarted  to  evaluate  the  re-acquisition  of  the  navigation 
signals  during  operational  conditions.  A  few  smaller  discontinuities  are  the 
result  of  satellite  switching  by  the  airborne  GPS  receiver.  During  satellite 
switching,  the  receiver  temporarily  enters  an  altitude-hold  mode.  The 
accuracy  of  the  resulting  differential  reference  trajectory  is  then  no  longer 
determined,  and  consequently,  the  trajectory  cannot  be  used  for  the  evaluation 
of  the  hybrid  GPS/LORAN  receiver.  Note  that  the  hybrid  receiver  can  still 
continue  to  provide  a  solution  based  on  the  three  remaining  satellites  and  two 
or  more  LORAN-C  transmitters. 

The  largest  position  errors  occur  during  the  middle  of  the  flight. 

These  deviations  are  caused  be  a  relatively  poor  GPS  geometry.  Also,  all 
sudden  changes  in  the  magnitude  of  the  two-dimensional  error  are  caused  by 
transitions  to  different  sets  of  four  GPS  satellites.  The  horizontal  position 
accuracy  for  the  hybrid  system,  based  on  all  measurements  (785  data  points), 
is  210  meters  (2  drms),  with  respect  to  the  DGPS  trajectory.  The  mean 
position  errors  in  the  North  and  East  directions  were  found  to  be  -52  meters 
and  30  meters,  respectively.  The  2  drms  positioning  accuracy  is  well  within 
all  current  requirements  for  enroute  navigation  (2778  m) ,  terminal  navigation 
(2037  m) ,  and  nonprecxsion  approaches  (556  m). 

Next,  the  performance  of  the  integrity  algorithm  was  addressed.  A 
straight-line  section  of  the  flight  data  was  used  which  lasted  for 
approximately  400  seconds.  During  this  time,  a  50  meters  per  second  ramp 
error  was  simulated  on  one  of  the  measurements,  as  depicted  in  Figure  9.3. 

The  actual  value  of  the  ramp  is  not  important  for  the  performance  of  the 
detection  algorithm,  since  the  algorithm  takes  "snapshots"  of  the  measurement 
data.  The  ramp  error  was  repeated  for  all  seven  measurements,  which  consisted 
of  four  GPS  satellites;  SV6,  SV9,  SVll,  and  SV12;  and  three  LORAN  stations: 
Carolina  Beach,  Dana,  and  Seneca.  Figures  9.4  and  9.5  show  the  resulting 
horizontal  radial  position  errors  and  the  integrity  or  range  residual 
parameters.  The  integrity  parameter  corresponds  to  e^  in  Figure  4.1. 

In  the  absence  of  malfunctions,  both  the  integrity  parameter  and  the 
horizontal  radial  position  error  are  well  below  200  m.  In  the  case  of  GPS  or 
Dana  malfunctions,  the  integrity  parameter  grows  faster  than  the  corresponding 
radial  position  error.  However,  both  Seneca  and  Carolina  Beach  malfunctions 
cause  the  radial  position  error  to  grow  rather  rapidly.  Detection  of  these 
malfunctions  is  still  possible  as  indicated  by  Figure  9.4, 
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Figure  9.1  Hybrid  GPS/LORAN  two-dimensional  position  errors. 
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Error  in  meters 


Figure  9.3  Sinulated  50  meters  per  second  ramp  failure. 
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Figure  9.4  Horizontal  radial  position  error  and  integrity  parameter 
for  a  simulated  50  meters  per  second  ramp  failure. 
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Figure  9.5  Horizontal  radial  position  error  and  integrity  parameter 
for  a  simulated  50  meters  per  second  ram.p  failure. 
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but  the  threshold  for  the  integrity  parameter  must  be  set  much  lower  than  for 
all  other  cases. 

Further  analysis  of  the  measurement  geometries  for  all  the  cases 
depicted  in  Figures  9.4  and  9.5  explains  the  larger  sensitivity  to  Seneca  and 
Carolina  Beach  malfunctions.  The  measurement  geometry  is  shown  in  Figure  9.6, 
where  the  arrows  represent  the  line-of-sight  vectors  to  the  GPS  satellites  and 
the  LORAN  transmitters  projected  onto  the  horizontal  plane  and  normalized. 

The  horizontal  dilutions  of  precision  (HDOP)  for  the  all-in-view  solution  and 
for  each  of  the  sub-sets  are  given  in  Table  9.1. 


I - 1 - 1 

I  Solution  I  HDOP  I 


I  All-in-view  |  1.08  | 
1  Sub-set  without  Seneca  |  1.94  j 
I  Sub-set  without  Carolina  Beach  |  2.17  | 
1  Sub-set  without  Dana  j  1.65  j 
1  Sub-set  without  Satellite  6  |  1.14  | 
I  Sub-set  without  Satellite  9  j  1.15  j 
j  Sub-set  without  Satellite  11  j  1.14  | 
j  Sub-set  without  Satellite  12  j  1.13  j 

I _ I _ I 


Table  9.1  Horizontal  dilutions  of  precision  (HDC?)  for  the  all-in-view 
solution  and  for  each  of  the  sub-seus. 

From  Figure  9.6  and  Table  9.1  it  can  be  seen  that  the  solutir>ns  without  Seneca 
and  Carolina  Beach  significantly  affect  the  remaining  geometry.  For  both 
cases,  the  HDOP  without  either  Seneca  or  Carolina  Beach  is  approximately  twice 
as  large  as  the  all-in-view  solution.  This  means  that  a  measurement  error  in 
either  of  these  signals  has  a  much  larger  effect  on  the  horizontal  radial 
position  error  than  any  of  the  other  five  signals,  which  is  clearly  shown  by 
Figures  9.4  and  9.5. 

In  addition  to  the  effects  of  geometry.  Figures  9.4  and  9.5  also  show 
the  effects  of  initial  bias  errors  in  the  measurements.  In  some  of  the  cases, 
the  integrity  parameter  initially  decreases,  which  indicates  that  the  injected 
failure  compensates  for  the  measurement  error.  As  expected,  this  also 
corresponds  to  a  smaller  radial  position  error.  After  10  to  20  seconds, 
however,  all  measurement  error  has  been  compensated  for  by  the  injected  error, 
and  the  integrity  parameter  will  steadily  grow. 
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10.  FLIGHT  TESTING  OF  THE  GPS/LORAN  RECEIVER 


The  prototype  hybrid  GPS/LORAN  receiver  has  been  flight  tested  on  two 
occasions  at  the  time  of  this  writing.  The  first  flight  test  occurred  on 
August  21,  1990;  and  the  second  flight  test  took  place  on  August  23,  1990. 

Both  flight  tests  were  performed  in  the  vicinity  of  the  Ohio  University 
airport  at  Albany,  Ohio.  The  duration  of  flight  number  one  was  approximately 
34  minutes,  and  the  duration  of  flight  number  two  was  approximately  52 
minutes.  Figures  10.1  and  10.2  show  the  ground  tracks  for  the  two  flight 
tests.  The  ground  tracks  reflect  the  simulated  signal  malfunctions;  both 
flight  tracks  contain  many  flight  path  deviations  caused  by  the  injection  of 
simulated  signal  malfunctions.  Tables  10.1  and  10.2  list  the  sequence  of 
events  for  each  of  the  flight  tests.  The  emphasis  of  the  flight  tests  was  on 
the  following  two  items: 

1.  Operational  verification  of  a  realtime  hybrid  GPS/LORAN  receiver. 

2.  Preliminary  assessment  of  the  flight  technical  error  and  the  impact 

of  failure  modes. 

Based  on  the  flight  tests,  it  was  apparent  that  the  hybrid  GPS/LORAN  receiver 
performed  in  accordance  with  its  design.  Both  test  pilots  noted  that  the 
course  deviation  indicator  is  very  responsive  and  that  the  indicated  course 
compares  favorably  with  other  area  navigation  equipment.  Also,  cross-checks 
with  DME,  VOR,  NDB,  and  the  ILS  localizer  indicated  that  the  GPS/LORAN 
horizontal  accuracy  is  on  the  order  of  0.1  nmi.  This  compares  well  with  the 
laboratory  results  provided  in  Chapter  8,  which  are  based  on  flight  data 
collected  under  similar  circumstances. 

The  flight  technical  error  is  the  accuracy  with  which  the  aircraft  is 
controlled  as  measured  by  the  indicated  aircraft  position  with  respect  to  the 
indicated  command  or  desired  position.  Only  the  cross-track  portion  of  the 
flight  technical  error  is  used  for  analysis  purposes,  (The  cross-track  error 
is  identical  to  the  signal  send  to  the  course  deviation  indicator).  The 
analysis  of  the  flight  technical  error  is  very  preliminary.  Because  of 
instrument  meteorological  conditions  (IMG)  during  take-off  and  landing,  the 
GPS/LORAN  course  could  not  be  flown.  Also,  both  pilots  were  briefed  about  the 
nature  of  the  flight  experiment  and  they  were  specifically  instructed  to  fly 
the  GDI  as  closely  as  possible.  Nevertheless,  the  FTE  error  traces  in  terms 
of  the  GDI  deviations  for  both  flights  are  shown  in  Figures  10.3  and  10.4. 

Note  that  for  flight  number  one,  the  full-scale  GDI  deflection  is  ±5  nmi 
(enroute  guidance)  or  ±1.25  nmi  (approach).  Both  test  pilots  found  a  smaller 
full-scale  GDI  deflection  nmi  desireable  for  the  enroute  mode.  Therefore, 
flight  number  two  used  ±2.5  nmi  for  the  enroute  guidance  and  ±1.25  for  the 
approach  mode.  Both  pilots  utilized  an  ILS  localizer  approach  to  runway  25  at 
the  end  of  the  mission.  As  a  result,  the  latter  part  of  the  FTE  traces  is  not 
representative  for  the  GPS/LORAN  FTE;  however,  it  is  a  good  indication  of  the 
GPS/LORAN  cross-track  error.  For  both  flights,  the  GPS/LORAN  cross-track 
error  is  of  the  order  of  0.1  nmi.  Since  the  offset  is  the  same  for  both 
flights,  it  is  most  likely  caused  by  uncorrected  LORAN  propagation  delays  and 
secondary  station  clock  offsets.  The  spikes  in  the  FTE  traces  are  either 
caused  by  waypoint  switching  or  by  the  injection  of  step  failures.  These  do 
not  reflect  the  ability  of  the  pilot  to  fly  the  indicated  GPS/LORAN  course. 
Because  of  the  capability  of  the  GPS  receiver  to  always  know  the  exact  value 
of  the  FTE,  a  study  is  recommended  to  Investigate  the  effects  of  flagging  a 
nonprecision  approach  based  on  a  FTE  of  such  magnitude  that  the  approach 
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Figure  10.1  Hybrid  GPS/LORAH  ground  track  for  flight  one,  August  21 
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gure  10.2  Hybrid  GPS/LORAH  ground  track  for  flight  two,  August  23,  1990. 
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Figure  10.3  Flight  technical  error  trace  for  flight  one. 
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Figure  10. A  Flight  technical  error  trace  for  flight  two. 


1  Time  (GMT)  |  Event 


17:09:00 

17:09:30 

17:10:37 

17:12:00 

17:U:A7 

17:15:00 

17:15:05 

17:15:57 

17:17:33 

17:17:45 


I  17:21:13 
i  17:24:35 

I  17:25:05 

i  17:27:35 
I  17:28:30 
I  17:29:00 
I  17:30:00 
I  17:30:10 
I  17:31:10 
I  17:31:23 
I  17:32:30 
I  17:34:15 
I  17:38:20 
I  17:38:50 
I  17:39:40 
I  17:41:15 

I  17:42:00 
I  17:43:00 
!  17:43:53 


Table  10. 


Take-off  from  UNI,  to  Henderson  VOR. 

GPS/LORAN  program  stopped  and  restarted  to  verify  correct 
operation  of  the  data  collection. 

Receiver  is  tracking  SV2,  SVll,  SV13,  SV14,  Seneca,  Dana, 
and  Carolina  Beach. 

GDI  switched  to  expanded  scale  per  pilot  request. 

Failure  1:  SV2  100  m/s  ramp. 

Integrity  threshold  of  300  m  exceeded,  GDI  is  flagged. 

Reset  failure  1. 

Failure  2:  SV2  25  m/s  ramp. 

Range  error  has  grown  to  2400  m,  integrity  parameter  =  316. 
Reset  failure  2. 

DME  cross-check,  Henderson  DME  agrees  to  within  0.1  nmi  with 
range  indicated  by  GPS/LORAN;  baro-altitude  -  4,000  ft. 
Standard  left  turn  towards  the  University  NDB. 

GPS/LORAN  program  stopped  and  restarted  to  save  collected 
data  to  flexible  disk. 

GPS/LORAN  receiver  operational,  waypoints:  from  Henderson 
VOR  to  University  NDB,  GDI  in  enroute  mode. 

Failure  3:  Carolina  Beach  10,000  m  step. 

Reset  failure  3. 

Failure  4:  SV14  50  m/s  ramp. 

Reset  failure  4. 

Failure  5:  Dana  100  m/s  ramp. 

Reset  failure  5. 

Initiation  of  descend. 

GDI  switched  to  expanded  scale  (approach  mode) . 

Waypoint  changed  to  threshold  of  RW  25. 

6.6  nmi  from  runway  threshold,  baro-altitude  -  2,500  ft 
Cross-check;  good  agreement  with  ILS  localizer. 

Passing  university  NDB,  4.76  nmi  from  threshold  of  RW  25. 
Cross-check  on  final;  GPS/LORAN  needle  shows  a  slight 
deflection  to  the  left. 

Visual  to  RW  25,  preparing  for  a  full  stop. 

Touch-down. 

GPS/LORAN  receiver  program  stopped. 


Sequence  of  events  for  flight  test  one,  August  21,  1990. 
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Time  (GMT) 


20:01:50  | 

20:0A:00  j 

i 

20:06:25  ] 

20:08:03  j 
20:08:20  1 
20:10:30  j 
20:11:07  j 
20:14:17  j 

20:16:00  | 
20:16:09  I 

20:20:00  | 

20:20:37  | 

20:22:23  | 

20:24:57  j 
20:27:25  | 

20:28:29  i 
20:28? ;0  I 

20:31:29  | 

20:32:25  i 
20:32:50  j 
20:35:10  j 
20:36:00  j 
20:38:07  | 

20:38:17  j 
20:40:01  j 
20:41:01  j 
20:43:43  j 
20:44:45  1 

20:46:30  | 

20:48:23  j 

20  :*49  :  11  1 

20:53:43  j 
20:54:43  I 


Event 


Take-off  from  UNI,  to  Yellow  Bud  VOR. 

At  2,500  ft  and  climbing.  Receiver  is  tracking  SV9,  SV13, 
SV14,  SV20,  Seneca,  Dana,  and  Carolina  Beach. 

Level  at  4,000  ft,  difficulty  tracking  SV20 
Failure  6:  SV9  10,000  m  step. 

Reset  failure  6. 

Failure  7:  Dana  25  m/s  ramp. 

GDI  is  flagged,  reset  of  failure  6. 

Failure  8:  Carolina  Beach  1,000  m  step,  little  effect. 
Difficulty  tracking  SV9.  Reset  failure  8. 

Failure  9:  SV20,  50  m/s  ramp. 

GDI  is  flagged,  still  difficulties  tracking  SV9 . 

Reset  failure  9. 

GPS/LORAN  program  stopped  and  restarted  to  save  collected 
data  to  flexible  disk. 

GPS/LORAN  receiver  operational,  waypoints  from  University 
NDB  to  Henderson  VOR.  Tracking  SV16,  SV13,  SV14,  SV20. 
Waypoints  from  Henderson  VOR  to  University  NDB. 

Passing  Yellow  Bud  VOR 

Failure  10:  SV14,  1,000  m  step.  GPS/LORAN  altitude  changed 
by  +3,000  ft. 

Reset  failure  10. 

Failure  11:  Seneca  10,000  m  step,  CDI  is  flagged, 
reset  failure  11. 

Failure  12:  SV16,  5  m/s  ramp. 

Reset  failure  12,  error  is  growing  too  slowly. 

Failure  13:  SV16,  10  m/s  ramp. 

CDI  is  flagged,  reset  failure  13,  CDI  changes  by  -30. 
Descend  to  3,000  ft  per  ATC  request,  16  nmi  to  beacon. 
Failure  14:  Dana,  10,000  m  step,  CDI  is  flagged. 

Reset  failure  14. 

Failure  15;  SV13,  100  m/s  ramp. 

CDI  is  flagged.  Reset  failure  15. 

8.3  nmi  to  NDB 

Switched  to  approach  mode,  waypoint  from  University  NDB  to 
threshold  of  RW  25. 

Procedure  turn  around  NDB. 

On  course.  Range  8.5  nmi  to  RW  25,  GPS/LORAN  CDI  needle 
slightly  to  the  left  compared  to  the  ILS  localizer. 

At  2,500  ft,  crossed  over  NDB. 

Touch-down. 

GPS/LORAN  receiver  program  stopped. 


Table  10.2  Sequence  of  events  for  flight  test  two,  August  23,  1990. 
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cannot  be  safely  completed.  Potential  benefits  are  increased  safety  and  a 
possible  reduction  in  the  size  of  obstacle  clearance  surfaces.  Furthermore, 
this  would  also  accomodate  the  development  of  generic  approach  procedures  for 
earth  referenced  navigation  systems. 

To  analyze  the  impact  of  signal  malfunctions,  a  total  of  fifteen 
simulated  failures  were  injected  into  the  measurement  data  during  the  two  test 
flights.  Both  sudden  errors  (step  failures)  and  slowly  increasing  errors 
(ramp  failures)  were  simulated.  Figures  10.5  and  10.6  show  the  horizontal 
radial  differences  between  the  unperturbed  GPS/LORAN  flight  trajectory  and  the 
trajectory  with  the  simulated  signal  failures,  as  well  as  the  integrity 
parameter  (range  residual  parameter).  Also  indicated  in  these  two  figures  are 
the  numbers  corresponding  to  each  of  the  failures.  Flight  number  one  used  an 
integrity  parameter  threshold  of  300  m,  flight  number  two  used  an  integrity 
parameter  threshold  of  400  m.  All  signal  malfunctions  which  would  have  caused 
unacceptable  course  deviations  were  detected  by  the  integrity  algorithm. 

Most  failures  increase  the  integrity  parameter  and  the  radial  position 
error  as  expected,  except  for  failure  number  8.  This  failure  is  a  1,000  m. 
step  on  the  Carolina  Beach  LORAN  measurement.  The  integrity  parameter  reduces 
from  300  meters  to  approximately  40  meters,  while  the  radial  position 
difference  grows  to  410  meters.  This  type  of  error  could  easily  be  detected 
by  a  time-history  filter  on  each  of  the  measurements.  Moreover,  the  Carolina 
Beach  transmitter  would  probably  have  "blinked"  under  the  simulated  condition. 
However,  if  this  error  would  be  slowly  increasing,  detection  would  be  much 
harder.  Closer  observation  of  the  data  showed  that  the  GPS  receiver 
experienced  difficulties  tracking  satellite  number  9.  As  a  result,  the 
Carolina  Beach  malfunction  is  expected  to  have  a  large  impact  on  the  position 
error  (see  also  Chapter  8).  Although  this  error  might  not  be  detected  at  the 
time  the  radial  error  threshold  is  exceeded,  the  least  squares  residual  method 
would  eventually  raise  the  flag. 

Although  the  test  pilots  knew  of  our  capability  to  inject  signal 
failures,  relevant  pilot  remarks  with  respect  to  each  of  the  fifteen  failure 
scenarios  are  given  in  Table  10.3.  It  is  interesting  to  note  that  both  test 
pilots  detected  significant  failures  before  the  integrity  threshold  was 
surpassed.  As  expected,  sudden  failures  are  easy  to  detect  by  the  pilot; 
however,  this  is  only  true  if  the  pilot  almost  continuously  monitors  the  CDI, 
otherwise,  small  step  failures  go  undetected  by  the  pilot. 

In  the  presence  of  strong  winds,  it  is  difficult  to  tell  the  difference 
between  a  heading  correction  to  compensate  for  a  cross  wind,  or  a  heading 
correction  caused  by  a  malfunctioning  navigation  signal,  especially  if  the  net 
effect  is  in  the  same  direction.  It  also  became  evident  that  the  pilots  were 
using  the  change  in  the  indicated  magnetic  course  to  detect  the  small  ramp 
failures.  In  addition  to  the  pilot  remarks,  the  GPS/LORAN  system  operators  in 
the  back  of  the  aircraft  also  noticed  large  climb /d'^scend  rates  and  large 
negative  altitudes  on  the  status  screen  as  a  result  of  the  injected  failures. 
In  some  of  the  cases,  the  unreasonableness  of  the  GPS/LORAN  altitude  was 
apparent  long  before  the  integrity  parameter  would  exceed  the  threshold. 

Based  on  the  results  of  the  flight  experiments,  the  following  is 
recommended  for  further  investigations:  evaluation  of  integrity  and  failure 
isolation  schemes  which  are  not  just  based  on  the  inconsistency  of  the 
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Figure  10.6  Horizontal  radial  difference  between  the  unperturbed  GPS/LORAN 
flight  trajectory  and  the  trajectory  with  simulated  signal  failures;  and  the 

integrity  parameter  for  flight  two. 
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Failure  Description 


Test  Pilot  Remark 


1.  SV2  100  m/s  ramp 


2.  SV2  25  m/s  ramp 


3.  Carolina  10,000  m  step 

4.  SV14  50  m/s  ramp 

5.  Dana  100  m/s  ramp 


After  a  few  seconds,  pilot  notices  a  deviation 
to  the  left,  pilot  reports  that  he  is 
correcting  for  the  wind.  Just  before  the  flag 
appears,  pilot  confirms  a  definite  command  to 
fly  left. 

Flag  is  noticed  immediately  by  pilot. 

After  approx.  30  seconds,  pilot  notices  a  slow 
deviation  to  the  left.  Flag  is  noticed  by  the 
pilot  within  one  second.  Pilot  reports  a 
heading  change  of  10  degrees. 

Pilot  notices  flag  immediately,  GDI  deviation 
is  not  significant. 

Error  discarded,  pilot  workload  did  not  allow 
for  monitoring/ flying  the  GPS/LORAN  GDI. 

Error  discarded  because  of  pilot  workload. 


6.  SV9  10,000  m  step 

7.  Dana  25  m/s  ramp 


8.  Carolina  1,000  m  step 

9.  SV20  50  m/s  ramp 

10.  SV14  1,000  m  step 

11.  Seneca  10,000  m  step 

12.  S'^16  5  m/s  ramp 

13.  SV16  10  m/s  ramp 


14.  Dana  10,000  m  step 

15.  SV13  100  m/s  ramp 


j  Pilot  notices  immediately  a  major  error. 

1  Error  not  noticed  by  pilot  until  several 
I  seconds  after  GDI  flag.  Pilot  remarks  that 
1  strong  winds  are  present,  he  thought  that  he 
1  was  correcting  for  a  cross  wind. 

I  Not  much  effect  on  the  navigation,  not 
I  noticed  by  pilot. 

I  GDI  flag  detected  by  pilot  within  a  few  sec. 
I  No  flag,  but  pilots  notices  wrong  GDI  track. 
I  GDI  flagged,  pilot  notices  error. 

I  Malfunction  not  noticed  by  pilot,  error 
I  injection  aborted,  it  would  take  too  long 
1  before  GDI  would  be  flagged, 
i  Pilot  starts  noticing  malfunction  when  the 
1  integrity  parameter  is  340  m.  GDI  flag  is 
I  detected  quickly. 

1  Error  discarded  because  of  pilot  workload. 

I  GDI  was  flagged. 

j  Pilot  notices  slight  jerk  to  the  right, 

I  which  he  notices  because  the  correction  is 
I  in  a  direction  opposite  the  wind. 


Table  10.3  GPS/LORAN  simulated  failure  scenarios  with  test  pilots’  remarks. 


GPS/LORAN  measurements,  but  which  also  take  the  following  information  into 
account : 


-  reasonableness  of  the  climb/descend  rate  and  altitude  as  indicated  by 
GPS/LORAN; 

-  rate  of  change  and  magnitude  of  the  differences  between  the  indicated 
GPS/LORAN  heading,  and  the  calculated  heading  and/or  indicated  magnetic 
heading. 
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1 1 .  CONCLUSIONS 


The  first  part  of  this  report  shows  the  derivation  of  the  ordinary  least 
squares  (OLS)  estimator  from  the  linear  model  using  the  Projection  Theorem. 

The  performance  of  the  OLS  is  subjected  to  a  detailed  error  analysis  which 
includes  the  effects  of  geometry  and  measurement  errors,  both  bias  and  noise. 
Using  similar  procedures,  the  extended  Ridge  and  Kalman  filters  are  derived 
which  are  representative  for  the  classes  of  biased  and  unbiased  estimators, 
respectively.  It  is  shown  that  in  the  presence  of  measurement  bias  errors, 
both  estimators  converge  to  the  OLS  inflation  of  the  position  error.  It  is 
therefore  concluded  that  in  the  case  of  dominant  measurement  bias  errors  (e.g. 
GPS  or  LORAN),  the  integrity  monitoring  performance  offered  by  the  OLS 
estimator  cannot  be  improved  upon.  It  was  also  found  that  the  Ridge  estimator 
can  be  used  to  explain/optimize  the  performance  of  a  mismatched  Kalman  filter. 
This  provides  a  very  helpful  insight  into  the  performance  of  the  Kalman  filter 
in  the  presence  of  unmodeled  dynamics. 

Based  on  the  conclusion  that  the  OLS  estimator  is  sufficient  to  perform 
integrity  monitoring,  a  general  solution  methodology  is  presented  for  a 
multisensor  navigation  system.  The  range  residual  integrity  parameter  was 
derived,  again  using  the  Projection  Theorem. 

The  above  results  have  been  applied  to  a  realtime  prototype  hybrid 
GPS /LORAN  receiver.  Two  flight  tests  revealed  that  the  hybrid  GPS /LORAN 
receiver  performs  in  accordance  with  its  design.  The  course  deviation 
indicator  is  responsive  and  the  indicated  course  compares  favorably  with  other 
area  navigation  equipment.  Out  of  a  total  of  fifteen  simulated  signal 
malfunctions,  twelve  malfunctions  which  would  have  caused  unacceptable  course 
deviations  were  detected  by  the  range  residual  integrity  algorithm. 

Equal  weighting  of  GPS  and  LORAN  measurements  is  used  for  the  flight 
tests.  Therefore,  the  accuracy  of  the  hybrid  system  will  be  mostly  determined 
by  the  LORAN  measurementc .  Standard  LORAN  propagation  models  are  used  such 
that  the  achieved  accuracies  are  representative  for  current  LORAN  receivers. 
Because  of  this,  the  accuracy  of  the  hybrid  system  will  not  be  as  good  as  that 
provided  by  GPS;  however,  the  availability  and  integrity  of  the  hybrid  system 
exceeds  that  of  GPS  by  several  orders  of  magnitude  [A].  At  the  same  time,  the 
hybrid  system  accuracies  were  still  found  to  be  well  within  all  current 
requirements  as  listed  in  Chapter  3.  It  is  also  concluded  that  the  accuracy 
of  hybrid  GPS/LORAN  can  be  improved  upon  significantly  through  LORAN 
calibration  or  by  using  a  weighting  matrix  in  the  hybrid  solution,  see  Section 
6.1. 


In  summary,  hybrid  GPS/LORAN  is  a  successful  case  study  of  fully  hybrid, 
multisensor  navigation.  Similar  performance  characteristics  may  be 
anticipated  for  other  multisensor  systems  based  on  sensors  such  as  GPS,  LORAN, 
GLONASS,  Omega,  and  baro-altimeter . 
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12.  RECOMMENDATIONS 


Based  on  the  findings  presented  in  this  report,  the  following 
recommendations  are  made  for  further  efforts  in  the  areas  of  multisensor 
navigation  and  GPS/LORAN: 

1 .  Development  of  a  generic  approach  procedure  for  earth  referenced 
navigation  systems. 

2.  Study  of  the  effects  of  flagging  a  nonprecision  approach  based  on  a 
flight  technical  error  of  such  magnitude  that  the  approach  cannot  be  safely 
completed.  Potential  benefits  are  Increased  safety  and  a  possible  reduction 
in  the  size  of  obstacle  clearance  surfaces. 

3.  Continued  integration  of  hybrid  GPS/LORAN  into  the  National  Airspace 
System.  The  results  of  this  study  can  also  be  applied  to  multisensor 
navigation  systems  in  general.  It  is  recommended  that  this  effort  takes  place 
along  the  following  avenues : 

-  Continued  flight  testing  of  the  prototype  hybrid  GPS/LORAN  receiver  to 
address  the  flight  technical  error  and  the  impact  of  failure  modes  on 
enroute  navigation  and  nonprecision  approach  operations; 

-  Detailed  study  of  failure  Isolation  schemes,  which  is  required  for 
sole  means  navigation; 

-  Development  and  evaluation  of  criteria  to  be  used  for  the 
certification  of  hybrid  GPS/LORAN  receivers; 

-  Development  and  evaluation  of  criteria  to  be  used  for  the  definition 
of  sole  means  of  navigation  systems. 

4.  Specific  efforts  to  further  integrate  GPS  and  LORAN  should  also 
address  promising  accuracy  improvement  techniques,  which  have  the  potential  to 
achieve  LORAN  measurement  accuracies  comparable  to  those  provided  by  GPS  [4, 

8,  61-63]; 

-  Use  of  a  weighting  matrix  W  to  incorporate  the  statistical  knowledge 
of  the  measurements  (see  Section  6.1); 

-  Calibration  of  LORAN  using  validated  GPS  positions; 

-  Use  of  improved  LORAN  propagation  models  which  could  contain  seasonal 
correction  data  based  on  the  LORAN-C  monitor  network. 

5.  Evaluation  of  integrity  and  failure  isolation  schemes  which  are  not 
just  based  on  the  inconsistency  of  the  GPS/LORAN  measurements,  but  which  also 
take  the  following  information  into  account: 

-  reasonableness  of  the  climb/descend  rate  and  altitude  as  Indicated  by 
GPS/LORAN; 

-  rate  of  change  and  magnitude  of  the  differences  between  the  Indicated 
GPS/LORAN  heading,  and  the  calculated  heading  and/or  indicated  magnetic 
heading. 
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